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Abstract 

This  document  serves  as  a  final  report  to  the  Air  Force  Office  of  Scientific  Research 
regarding  activities  performed  under  Grant  #  FA9550-09-1-0069,  “Advanced  Chemical 
Modeling  For  Turbulent  Combustion  Simulations.”  The  project  had  two  goals.  The 
first  was  the  development  of  new  sub-filter  models  for  large  eddy  simulation  of  tur¬ 
bulent  combustion,  and  the  second  was  the  development  of  chemical  mechanisms  for 
jet  fuel  surrogates.  This  report  presents  and  discusses  the  models  and  the  results  that 
were  produced  within  each  project  component.  The  sub-filter  modeling  work  partic¬ 
ularly  focuses  on  the  development  of  a  framework  for  describing  multiple  combustion 
regimes  using  the  flamelet  approach,  on  describing  the  scalar  dissipation  rate  in  tur¬ 
bulent  non-premixed  combustion,  and  on  modeling  strain  effects  in  turbulent  premixed 
combustion.  The  chemical  mechanism  work  proposes  a  method  for  defining  jet  fuel 
surrogates,  describes  how  sub-mechanisms  for  a  variety  of  alkanes  and  aromatics  can 
be  incorporated  into  the  definition,  and  finally  creates  and  validates  a  mechanism  that 
serves  as  an  accurate  surrogate  for  jet  fuel  behavior. 
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1  Introduction 


The  application  of  computational  models  to  multi-physics  flow  problems  is  becoming  in¬ 
creasingly  routine.  The  proliferation  of  computational  resources,  flow  physics  software, 
and  data  processing  capabilities  has  provided  modelers  with  the  ability  to  analyze  realis¬ 
tic  engineering  systems  that  are  characterized  by  complex  geometry,  multiphase  flow,  and 
non-linear  chemical  interaction.  In  spite  of  their  availability,  computational  flow  physics 
models  produce  results  that  remain  subject  to  large  uncertainties.  These  uncertainties  neg¬ 
atively  impact  the  extent  to  which  the  model  results  can  be  relied  upon  for  engineering 
development  and  analysis.  A  variety  of  research  efforts  are  directed  toward  reducing  this 
uncertainty  and  improving  the  usefulness  of  flow  physics  modeling  frameworks.  Efforts  to 
develop  more  physically  accurate  models  of  reactive  flow  processes  are  among  the  most 
important  of  these  efforts.  Without  the  development  of  increasingly  accurate  physical  mod¬ 
els,  the  confidence  bounds  associated  with  multi-physics  engineering  analyses  cannot  be 
decreased. 

In  this  project  two  categories  of  physical  models  are  targeted.  The  first  category  consists 
of  descriptions  of  turbulence  and  chemistry  sub-filter  interactions  in  Large  Eddy  Simulation 
(LES).  The  second  category  consists  of  descriptions  of  the  evolution  of  the  chemical  kinetics 
of  jet  fuel.  While  both  of  these  modeling  topics  have  been  widely  studied  [1,  2,  3],  they 
represent  two  of  the  greatest  modeling  challenges  within  the  wider  landscape  of  physical  flow 
modeling.  For  example,  making  assumptions  about  the  combustion  regimes  that  are  found 
on  sub-filter  LES  scales  may  result  in  dramatic  mis- predictions  of  CO  [4],  while  uncertainties 
in  the  prescription  of  kinetics  rates  can  lead  to  vastly  different  NO  predictions  [5].  Advancing 
and  further  developing  these  models  requires  the  application  of  fundamental  theory,  physical 
insight,  and  detailed  comparisons  with  numerical  and  experimental  validation  cases.  Each 
of  these  approaches  will  be  relied  upon  in  the  sections  that  follow. 

Activities  related  to  the  first  project  topic,  which  is  the  development  and  validation  of  a 
series  of  sub- filter  models  for  reactive  LES,  are  discussed  in  section  2.  This  section  addresses 
modeling  challenges  from  across  the  entire  spectrum  of  combustion  physics  that  are  encoun¬ 
tered  in  reacting  flows.  The  greatest  of  these  challenges  is  the  treatment  of  the  assumptions 
that  are  inherent  in  the  flamelet  approach.  This  approach  is  used  for  combustion  modeling 
because  it  has  historically  offered  the  most  accurate  means  of  describing  unresolved  turbu¬ 
lence  and  chemistry  interactions  in  LES.  Flamelets,  however,  make  detailed  assumptions 
about  the  physics  that  govern  combustion.  These  assumptions  must  change  when  the  lo¬ 
cal  flow  conditions  change  and  give  rise  to  different  types  of  combustion  processes.  Each 
combustion  process  can  be  characterized  as  belonging  to  a  particular  combustion  regime, 
and  here  a  model  is  developed  that  uses  local  transport  information  from  a  flow  solver  to 
predict  which  combustion  regime  will  be  present  at  any  point  in  a  simulation  domain.  Once 
this  regime  is  known,  appropriate  flamelet  assumptions  can  be  applied  in  the  flow  solver. 
Application  of  these  assumptions,  however,  requires  information  about  the  local  flow  held. 
In  the  non-premixed  combustion  regime,  the  most  critical  parameter  describing  the  local 
how  held  is  the  scalar  dissipation  rate.  While  widely  studied,  it  is  demonstrated  in  sec- 
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tion  2  that  this  quantity  continues  to  be  poorly  predicted  by  current  models.  In  response, 
a  new  model  that  adopts  traditional  LES  dynamic  machinery  in  a  novel  fashion  is  devel¬ 
oped  for  the  dissipation  rate.  Finally,  the  limit  of  premixed  combustion,  which  represents 
a  second  distinct  regime,  is  considered.  A  variety  of  models  describing  turbulent  premixed 
flame  propagation  are  available  in  literature.  The  most  widely  used  models,  however,  are 
characterized  by  an  inability  to  describe  flame  physics  when  turbulence  becomes  so  intense 
that  it  perturbs  the  flame’s  structure.  These  scenarios  (high  Karlovitz  number  flames)  are 
encountered  just  before  a  premixed  flame  is  extinguished  by  turbulence.  Describing  the 
underlying  physics  correctly  is  therefore  a  requirement  of  any  model  that  seeks  to  predict 
next  generation,  low  pollutant  premixed  engine  designs.  Here,  this  challenge  is  addressed 
by  formulating  a  new  method  of  describing  flame  structure  perturbations  in  LES. 

The  objective  of  the  second  project  topic  is  to  develop  a  capability  to  reliably  predict 
the  combustion  characteristics  of  fuel  oxidation  and  pollutant  emissions  from  engines.  The 
relevant  fuel  chemistry  must  be  accurately  modeled  if  this  goal  is  to  be  achieved.  Jet  fuel, 
however,  is  made  up  of  hundreds  of  hydrocarbons  in  proportions  that  can  vary  significantly 
between  fuel  samples.  Consequently,  only  average  fuel  properties  are  known  under  the  best 
of  circumstances.  A  simplified  fuel  representation  is  required  to  circumvent  this  problem. 
Typically,  jet  fuels  are  modeled  using  a  representative  surrogate  mixture,  which  is  simply  a 
well-defined  mixture  comprised  of  a  few  components.  This  surrogate’s  composition  is  chosen 
so  that  it  mimics  the  physical  and  chemical  properties  of  the  real  fuel  under  consideration. 
Guidelines  have  been  developed  by  surrogate  working  groups  [6,  7,  8]  to  determine  the 
individual  components  that  could  make  up  such  surrogate  mixtures  for  a  variety  of  specific 
fuels.  Given  knowledge  of  the  target  properties  of  the  real  fuel,  and  given  a  particular 
choice  of  individual  components  that  will  be  used  in  the  surrogate  mixture,  a  constrained 
optimization  approach  can  be  applied  to  determine  the  exact  composition  of  the  surrogate 
mixture.  An  optimization  tool  that  performs  this  task  was  developed  under  an  earlier 
AFOSR  grant,  and  is  used  in  this  project  to  determine  a  surrogate  mixture  composition 
that  best  matches  the  targeted  fuel  properties. 

Once  this  optimization  procedure  was  developed,  the  critical  modeling  need  shifted  to¬ 
ward  the  design  of  compact  and  reliable  kinetic  models  that  accurately  describe  the  behavior 
of  the  components  used  in  defining  the  surrogate.  This  issue  is  addressed  in  section  3.  The 
kinetic  models  that  are  needed  must  be  comprehensively  validated  for  the  individual  compo¬ 
nents  of  the  surrogate  fuel,  and  then  again  validated  after  they  are  combined  to  form  a  final 
surrogate  mixture.  Though  detailed  mechanisms  for  many  fuel  components  are  available 
in  literature,  these  existing  mechanisms  typically  consist  of  several  thousand  species  and 
reactions.  Even  single  component  mechanisms  can  reach  this  level  of  complexity,  which  pre¬ 
cludes  the  coupling  of  the  mechanism  with  a  tractable  flow  modeling  approach.  Moreover, 
component  mechanisms  from  different  sources  typically  describe  reaction  pathways  using 
different  methodologies.  These  disparate  descriptions  are  attributable  to  the  uncertainties 
that  exist  in  kinetic  rate  data.  Multiply  defined  pathways  render  the  process  of  combin¬ 
ing  detailed  component  mechanisms  nearly  impossible.  In  response  to  these  challenges,  a 
framework  has  been  developed  that  consists  of  first  reducing  individual  component  mecha- 
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nisms,  then  combining  them  into  a  mixture  characterizing  a  jet  fuel  surrogate,  and  finally 
further  reducing  and  validating  the  combined  surrogate  definition.  During  this  process, 
consistency  is  enforced  amongst  the  formulations  of  each  component  mechanism,  so  that 
the  kinetic  rates  of  each  reaction  class  do  not  suffer  from  duplication  or  conflicts  in  the  final 
mechanism. 
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2 


Advanced  Sub-Filter  Models  For  Large  Eddy  Simulation 
Of  Turbulent  Combustion 


2.1  Review  Of  Objectives 

Part  one  of  this  project  sought  to  develop  advanced  sub-filter  models  for  Large  Eddy  Sim¬ 
ulation  (LES)  of  turbulent  combustion.  Four  specific  goals  were  identified.  The  first  was 
the  development  of  a  flamelet  based  chemical  source  term  description  that  could  accurately 
predict  both  premixed  and  non-premixed  combustion  regimes.  The  second  goal  was  the 
creation  of  an  LES  description  of  premixed  flame  structure  that,  when  used  in  multi-regime 
combustion  environments,  could  seamlessly  integrate  with  typical  non-premixed  combus¬ 
tion  models.  Increasing  the  fidelity  of  flamelet  based  LES  predictions  of  pollutants  such  as 
NOx.  and  CO  was  a  third  goal.  Finally,  the  fourth  goal  was  the  development  of  improved 
descriptions  of  turbulent  premixed  flame  propagation  speeds. 

Here,  the  modeling  and  theory  that  was  developed  to  accomplish  these  goals  is  high¬ 
lighted  in  detail.  First,  a  multi-regime  combustion  model  that  provides  for  the  specifica¬ 
tion  of  chemical  source  terms  in  different  regimes  is  dicsussed.  This  model  builds  upon  the 
group's  prior  regime  identification  work  to  provide  a  model  implementation  characterized  by 
improved  tractability.  The  model  is  reviewed  in  section  2.2.  Next,  this  multi-regime  model 
is  validated  using  a  canonical  flame  case:  a  partially  premixed  triple  flame.  These  cases  are 
described  and  analyzed  in  section  2.3.  They  demonstrate  that  the  multi-regime  approach  is 
capable  of  predicting  flame  behavior  in  both  premixed  and  non-premixed  regimes  without 
the  need  for  a  user  to  make  a  priori  assumptions  about  mixing  dynamics.  Additionally,  the 
validation  work  in  section  2.3  demonstrates  the  sensitivity  that  pollutant  predictions  have 
to  burning  regime.  This  sensitivity  implies  that  the  goal  of  increasing  pollutant  prediction 
accuracy  is  being  implicitly  accomplished  through  the  multi-regime  model. 

Most  combustion  modeling  approaches  require  accurate  predictions  of  the  sub-filter 
scalar  dissipation  rate.  This  quantity  directly  affects  the  evolution  of  pollutants,  flame 
stability,  and  heat  release.  In  section  2.4,  large  eddy  simulations  of  a  high  fidelity  direct 
numerical  simulation  (DNS)  flame  database  demonstrate  that  the  most  traditional  LES 
dissipation  rate  model  can  lead  to  dramatic  underpredictions  of  sub- filter  scalar  dissipation. 
These  underpredictions  are  addressed  by  considering  a  new  transport  equation  closure  of  the 
dissipation  rate.  Specifically,  an  additional  DNS  describing  scalar  mixing  in  homogeneous 
turbulence  is  used  to  analyze  closure  models  for  the  relevant  transport  equation.  This 
analysis  leads  to  the  development  of  a  new  adapted  dynamic  procedure.  As  shown  at 
the  end  of  section  2.4,  the  use  of  the  new  transport  equation  closure  in  LES  is  shown 
to  significantly  improve  dissipation  rate  predictions  relative  to  the  traditional  modeling 
approach.  As  a  result,  the  lift-off  height  and  pollutant  statistics  of  the  high  fidelity  flame 
DNS  are  reproduced  in  LES  with  much  greater  accuracy. 

Finally,  the  fourth  goal  of  turbulent  flame  speed  modeling  is  addressed  in  section  2.5.  A 
focus  is  placed  on  described  premixed  flame  propagation  in  high  Karlovitz  number  regimes, 
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where  turbulence  can  penetrate  the  inner  layer  of  a  flame  and  affect  its  structure.  This 
physical  process  has  been  historically  difficult  to  describe  in  LES  because  it  occurs  entirely 
on  unresolved  scales.  Traditional  flamelet  models  that  assume  a  steady  premixed  flame 
structure  fail  to  predict  flame  propagation  under  these  circumstances.  This  challenge  is 
addressed  by  introducing  a  new  tabulation  coordinate  in  the  flamelet  formalism  to  describe 
the  sub-filter  effects  that  turbulence  has  on  the  flame  front.  This  new  coordinate  is  shown  to 
map  out  a  wide  range  of  premixed  flame  structures.  The  resulting  flamelet  model  is  tested 
by  performing  LES  of  an  additional  high  fidelity  premixed  DNS  case.  The  new  model  is  able 
to  reproduce  the  statistics  and  structure  of  the  DNS  flame  where  the  traditional  premixed 
flamelet  approach  could  not. 

2.2  Regime  Independent  Sub-Filter  Combustion  Modeling 

Flamelet  combustion  models  for  LES  are  most  often  formulated  for  single  regimes.  Advanc¬ 
ing  the  predictive  capabilities  of  these  models  requires  the  development  of  formulations  that 
can  determine  which  mode  of  combustion  occurs  at  a  local  point  within  a  flow  solver.  This 
information  can  then  be  used  to  map  in  flamelet  solutions  whose  assumptions  are  consis¬ 
tent  with  the  local  combustion  regime.  In  this  section,  a  tractable  implementation  of  such 
a  multi-regime  flamelet  model  is  developed.  Such  a  discussion  must  begin  with  a  review  of 
single  regime  flamelet  implementations. 


2.2.1  Single  Regime  Flamelet  Models 

Burning  in  the  purely  non-premixed  regime  can  be  described  using  a  standard  non-premixed 
flamelet  model  [1,  9,  10,  11].  This  approach  relies  on  solutions  of  the  steady  non-premixed 
flamelet  equations,  which  can  be  written  [12] 

Xz  4>k  d2W  xz  ( x  i  (t>k(l>md2W\  xzd2(j)k  .  m 

■hf  w+'yl  *8?*  w )-pTdp  = 


X^^dcpdT  xz  v  (d$m  cj>mdw\  f  cp,m\  dT  xz  o2t 

P  2  Cp  dz  dz  +  P  2  ^  V  dZ  +  W  dZ  )  \  cp  J  dZ  P  2  dZ2 


*N.(2) 

Cp 


The  scalar  dissipation  rate  in  Eq.  (2),  xz  =  2 T>z(VZ  -VZ),  is  an  external  parameter  in  the 
equation  and  is  modeled  as  a  prescribed  function  of  Z  [1,  12], 


exp  ^—2  [erfc  x(  2 Z  )]2^ 
exp  (-2  [erfc-1(2Zst)]2) 


(3) 


where  the  reference  dissipation  rate  at  the  stoichiometric  mixture  fraction  is  xz,st-  The 
flamelet  space  boundary  conditions  that  are  relevant  for  the  laminar  flames  being  considered 
below  are  T=  300  K  at  both  Z= 0  and  Z=l. 

The  non-premixed  flamelet  solutions  can  be  tabulated  in  a  database  as  a  function  of  Z 
and  C  [10].  Representative  non-premixed  flamelets  that  are  used  for  modeling  the  validation 


cases  in  section  2.3,  and  their  associated  S-shaped  curve  [1],  are  plotted  in  Fig.  1.  The  plot  of 
the  S-shaped  curve  shows  that  several  flamelets  from  the  so-called  lower  branch  are  included 
in  the  flamelet  solution  set.  The  lower  branch  is  defined  with  respect  to  the  turning  point, 
or  the  maximum  reference  scalar  dissipation  rate  that  can  be  observed  in  a  flamelet  before 
quenching  occurs.  Flamelets  that  have  maximum  temperatures  and  dissipation  rates  that 
are  smaller  than  the  temperature  and  dissipation  rate  of  the  flamelet  at  the  turning  point 
are  said  to  exist  on  the  lower  branch  of  the  S-shaped  curve.  These  solutions  are  unstable  in 
the  sense  that  they  will  respond  to  small  perturbations  by  moving  either  toward  the  stable 
upper  branch  or  toward  a  stable  quenched  solution.  Because  of  this  response,  the  usefulness 
of  lower  branch  solutions  in  non-premixed  flamelet  modeling  remains  ambiguous. 

Compounding  this  ambiguity,  it  is  known  [13]  that  the  treatment  of  flamelet  space  below 
the  turning  point  can  influence  flame  behavior.  For  example,  when  the  flamelet  solution  at 
the  turning  point  is  taken  to  be  the  first  available  burning  solution  above  the  mixing  line 
(the  completely  quenched  state),  then  all  flamelet  information  between  the  mixing  line  and 
this  first  available  flamelet  would  be  determined  using  interpolation.  These  interpolated 
values  are  expected  to  be  different  than  the  values  associated  with  flamelet  solutions  from 
the  lower  branch  of  the  S-shaped  curve.  In  comparison  with  a  purely  non-premixed  flamelet 
model,  one  advantage  of  a  multi-regime  modeling  approach  is  that  it  tends  to  access  the 
premixed  regime  in  these  regions  of  Z  and  C  space  where  the  progress  variable  is  not 
near  its  maximum  value.  A  multi-regime  approach  therefore  negates  some  of  the  ambiguity 
associated  with  the  lower  branch  solutions.  Note  that  when  the  limit  of  purely  non-premixed 
combustion  is  considered,  the  lower  branch  flamelets  will  be  used  without  modification.  The 
fully  quenched  mixing  line  will  also  be  included  in  flamelet  space,  and  regions  between  the 
lower  branch  and  the  mixing  line  will  be  accessed  using  interpolation.  Once  non-premixed 
flamelet  databases  are  formed,  the  purely  non-premixed  model  can  be  applied  to  laminar 
flames  by  solving  transport  equations  for  the  mixture  fraction  Z  and  the  progress  variable 
C. 

A  purely  premixed  flamelet  model  can  be  similarly  developed.  The  premixed  model  is 
based  on  a  coupled  level  set  and  progress  variable  approach  [14]  that  is  extended  through 
the  use  of  a  technique  from  the  Flamelet  Generated  Manifold  (FGM)  methodology  [15]. 
Flamelet  solutions  are  created  for  this  purely  premixed  approach  by  using  the  FlameMaster 
program  [16]  to  solve  the  steady  unstrained  premixed  flamelet  equations, 
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where  s l,u  is  the  eigenvalue  of  the  problem  that  represents  the  laminar  flame  speed,  pu  is 
the  density  on  the  unburned  side  of  the  flame,  and  x  is  a  physical  space  coordinate.  The 
diffusion  velocity  V^%P  is  usually  modeled  by  assuming  Fickian  diffusion.  Representative 
premixed  flamelet  solutions  that  are  used  for  validation  in  section  2.3  are  plotted  in  Fig.  2. 
The  lean  and  rich  flammability  limits  associated  with  premixed  burning  appear  in  this 
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Figure  1:  Representative  on-premixed  flamelets.  First  plot:  flamelet  temperature  profiles;  Second  plot: 
flamelet  CO  mass  fraction  profiles;  Third  plot:  maximum  flamelet  temperature  vs.  \z, st;  Fourth  plot: 
maximum  NO  production  vs.  \z, st-  The  Xz,st  values  in  the  lower  row  plots  are  taken  from  the  flamelets  in 
the  upper  row. 


figure  at  mixture  fractions  of  approximately  Zlfl  =  0.03  and  Zrfl  =  0.26,  respectively. 
A  comparison  of  Figs.  1  and  2  reveals  the  effect  of  a  change  of  the  combustion  regime. 
For  example,  CO  mass  fractions  in  the  premixed  regime  reach  values  approximately  twice 
as  large  as  the  CO  mass  fractions  in  the  non-premixed  regime.  Similarly,  the  maximum 
observed  temperature  in  the  two  regimes  at  Z= 0.2  differs  by  approximately  150  K. 

Just  as  in  the  non-premixed  model,  the  flamelets  are  tabulated  as  a  function  of  Z 
and  C.  Any  scalar  quantity  4>k  can  then  be  accessed  as  4>k  =  4>k(Z,C).  Because  of  the 
premixed  flammability  limits,  however,  mixture  fraction  values  outside  of  Zlfl  and  Zrfl 
are  not  available  in  the  premixed  solution  space.  A  special  extension  procedure  is  therefore 
needed  to  describe  very  lean  or  very  rich  mixture  fractions.  Here,  a  technique  based  on 
the  flamelet  generated  manifolds  (FGM)  model  [15]  will  be  used  to  perform  the  extension. 
In  this  technique,  scalar  information  at  a  Z  value  sitting  between  a  flammability  limit  and 
a  Z  boundary  value  is  determined  by  interpolating  between  the  premixed  solution  at  the 
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Figure  2:  Representative  premixed  flamelets.  First  plot:  flamelet  temperature  profiles;  Second  plot:  flamelet 
CO  mass  fraction  profiles;  Third  plot:  burning  velocity  as  a  function  of  Z\  Fourth  plot:  flamelet  NO 
production  profiles. 


flammability  limit  and  the  unburned  solution  at  the  Z  boundary, 

MZ,C)  =  (, hiZRFL,C )  •  Zrft^J 

+  4>k(z  =  i, c  =  o)  •  ^ ^  v  z>  zRFL, 

Mz,  c)  =  MzLFl,c)-^z^~°_^ 

+  mz  =  o, c  =  o)  •  (zLlffllZq)  v  z<  Zlfl' 


(6) 

(7) 


Note  that  a  multi-regime  model  will  again  help  to  remove  the  ambiguity  associated  with 
this  extrapolation  because  it  tends  to  access  non-premixed  flamelet  space  at  mixture  frac¬ 
tion  values  outside  of  the  flammability  limits.  In  multi-regime  approaches,  the  particular 
interpolation  technique  that  is  used  for  premixed  flamelet  extension  takes  on  a  reduced 
importance. 

In  fully  resolved  simulations  such  as  the  laminar  flames  considered  here,  a  progress 
variable  equation  can  be  solved  without  the  use  of  any  special  numerical  techniques.  For 
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example,  a  progress  variable  is  solved  for  directly  in  the  counter  flow  flame  simulation 
presented  in  section  ??.  In  this  approach  a  progress  variable  source  term  from  a  flamelet 
solution  is  directly  introduced  in  the  G  transport  equation.  In  standard  implicitly  filtered 
LES  simulations,  however,  progress  variable  profiles  tend  to  be  underresolved  around  flame 
fronts,  and  special  numerical  techniques  are  needed  to  address  the  resulting  problems  [14, 
17,  18].  Therefore,  in  an  effort  to  consider  flamelet  combustion  models  that  are  applicable 
to  LES,  the  triple  flame  simulation  will  couple  the  solution  of  a  level  set  to  the  progress 
variable  held.  Although  it  is  not  needed  in  a  fully  resolved  simulation,  this  coupling  would 
ensure  that  the  progress  variable’s  description  of  the  flame  front  would  propagate  at  the 
correct  speed  even  if  the  progress  variable  held  were  to  transition  from  unburned  to  burned 
over  a  single  mesh  cell.  A  level  set  will  be  defined  using  a  held  variable  labeled  as  G,  and 
the  particular  level  set  of  this  variable  that  denotes  the  flame  front  is  set  as  G  =  Go .  This 
level  set  is  governed  by  the  equation  [1,  14,  17], 

r\  r\ 

-(G) +  %  —  (G)  =  ^(Vuk  +  sl,u)\VG\  V  G  =  G0, 

|VG|  =  1  V  G,  (8) 

where  k  is  the  curvature  of  the  hame  front,  where  V  is  the  progress  variable  diffusivity 
at  the  hame  front,  and  where  the  u  subscript  represents  conditioning  on  the  unburnt  side 
of  the  front.  When  solving  Eq.  (8),  the  hame  speed  sl,u  is  determined  by  accessing  the 
tabulated  premixed  flamelet  solutions.  The  computed  burning  velocities  shown  in  Fig.  2 
are  therefore  directly  incorporated  in  the  premixed  model  framework  through  the  level  set 
equation. 

In  the  premixed  hamelet  model,  the  transport  equation  for  the  level  set  is  solved  in 
addition  to  the  Z  and  G  expressions  shown  in  Eqs.  (??)  -  (??).  The  progress  variable 
equation  is  coupled  to  the  level  set  by  formulating  the  progress  variable  source  term,  ujc, 
as  a  function  of  both  the  level  set  and  the  tabulated  premixed  solution.  This  coupling 
ensures  that  the  iso-contour  of  the  progress  variable  that  is  associated  with  the  flame  front 
is  always  consistent  with  the  position  of  the  level  set.  The  details  of  the  coupling  are 
described  elsewhere  [14]. 

2.2.2  Multi-Regime  Flamelet  Model 

The  two  single  regime  models  that  have  been  presented  can  be  combined  to  form  a  multi¬ 
regime  flamelet  model.  This  multi-regime  model  must  accomplish  two  tasks.  It  must  first 
provide  a  means  of  distinguishing  which  regime  is  present  at  a  particular  mesh  cell  and 
particular  point  in  time,  and  it  must  additionally  provide  a  means  of  transitioning  between 
solutions  from  the  individual  regimes.  The  multi-regime  approach  that  is  used  here  is  based 
on  a  flamelet  transformation  in  which  the  statistical  dependence  of  mixture  fraction  and 
progress  variable  are  accounted  for  [19].  This  approach  will  now  be  described. 

In  the  asymptotic  limit  of  purely  non-premixed  combustion,  reactive  source  terms  are 
balanced  by  diffusive  transport  that  occurs  in  the  direction  of  mixture  fraction  gradients. 
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The  non-premixed  flamelet  equations  that  are  shown  in  Eqs.  (1)  -  (2)  are  derived  by  per¬ 
forming  a  coordinate  transformation  that  explicitly  captures  this  balance.  Conversely,  in 
the  limit  of  purely  premixed  combustion  reactive  source  terms  are  balanced  by  diffusive 
transport  that  occurs  along  vectors  where  the  mixture  fraction  is  constant.  The  premixed 
flamelet  equations  shown  in  Eq.  (5)  explicitly  capture  this  second  type  of  reaction  and  dif¬ 
fusion  balance.  If  a  multi-regime  model  is  to  determine  which  of  these  asymptotic  regimes 
is  most  descriptive  of  the  combustion  at  a  given  flow  location,  the  balance  between  reaction 
and  diffusion  processes  must  be  examined.  If  reaction  is  primarily  balanced  by  transport  of 
the  mixture  fraction,  then  the  local  regime  can  be  said  to  be  non-premixed.  Conversely,  if 
reaction  is  balanced  by  transport  along  vectors  of  constant  mixture  fraction,  then  the  local 
regime  can  be  said  to  be  premixed. 

Traditional  mixture  fraction  and  progress  variable  based  flamelet  methods  present  a 
challenge  to  the  examination  of  these  transport  processes.  This  challenge  arises  because 
Z  and  C  are  not  statistically  independent.  In  the  limit  of  non-premixed  combustion,  for 
example,  a  change  in  the  mixture  fraction  field  will  be  accompanied  by  a  change  in  the 
progress  variable  field.  Consequently,  observations  of  diffusive  progress  variable  transport 
do  not  conclusively  indicate  the  presence  of  a  particular  combustion  regime.  If  diffusive 
transport  is  to  be  definitively  associated  with  the  premixed  or  the  non-premixed  regime, 
the  statistical  dependence  of  C  and  Z  must  be  properly  accounted  for.  In  reference  [19], 
this  dependence  is  treated  by  defining  a  variable,  A,  that  is  analogous  to  the  traditional 
progress  variable.  The  A  variable  has  been  previously  considered  [10,  11],  and  can  be  most 
easily  described  as  a  modified  reactive  coordinate  whose  value  is  constant  over  a  single  non- 
premixed  flamelet.  Mathematically,  the  A  value  that  is  associated  with  a  single  flamelet 
might  be  defined  as  the  value  of  the  progress  variable  C  that  occurs  at  the  flamelet’s 
stoichiometric  mixture  fraction  [19].  Because  non-premixed  flamelet  solutions  each  have 
slightly  different  progress  variable  profiles,  the  A  value  associated  with  each  non-premixed 
flamelet  will  be  unique.  Additionally,  because  A  is  constant  over  an  entire  non-premixed 
flamelet,  it  provides  a  measure  of  reaction  progress  that  is  statistically  independent  of  the 
mixture  fraction.  A  is  therefore  useful  to  an  examination  of  the  nature  of  local  transport 
processes.  The  purely  premixed  limit  will  be  indicated  by  significant  transport  of  the  A 
variable  and  no  transport  of  the  Z  variable.  Conversely,  the  purely  non-premixed  limit 
will  be  indicated  by  significant  transport  of  the  Z  variable  and,  since  A  is  constant  along  a 
non-premixed  flamelet,  no  transport  of  the  A  variable. 

Figure  3  shows  a  series  of  non-premixed  flamelet  solutions  in  a  way  that  emphasizes  the 
physical  meaning  of  the  A  variable.  The  same  steady,  non-premixed  flamelet  solutions  that 
were  shown  in  Fig.  1  are  plotted  in  Fig.  3  in  (Z, A)  space.  The  height  of  the  plotted  surface 
indicates  the  local  value  of  the  progress  variable  C  in  this  space,  while  the  color  contours 
indicate  the  corresponding  value  of  temperature.  As  described  above,  the  A  variable  is 
constant  in  value  throughout  a  single  non-premixed  flamelet  solution.  Consequently,  a 
single  flamelet  solution  can  be  recovered  from  Fig.  3  by  extracting  data  along  a  constant 
value  of  the  A  coordinate.  Additionally,  Fig.  3  demonstrates  that  it  is  possible  to  move 
within  a  flamelet  solution  where  A  is  constant  and  still  allow  both  Z  and  C  to  change. 
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This  coordinate  independence  is  what  enables  the  rigorous  separation  of  non-premixed  and 
premixed  transport  processes.  One  final  aspect  of  Fig.  3  that  is  relevant  to  the  proposed 
multi-regime  modeling  approach  is  the  sharp  change  in  topology  that  occurs  along  the  line 
where  A  =  0.15.  This  topological  change  appears  as  a  discontinuity  in  the  slope  of  the 
progress  variable  surface,  and  it  identifies  the  last  steady,  burning  flamelet  solution  that 
is  solved  for.  The  only  flamelet  solution  with  a  lower  maximum  temperature  is  the  non¬ 
burning  mixing  solution  that  exists  at  A  =  0.0.  Fig.  3  demonstrates  how  linear  interpolation 
is  used  to  populate  this  unstable  region  of  flamelet  space. 
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Figure  3:  The  steady  non-premixed  flamelet  solutions  from  Fig.  1,  plotted  in  (Z, A)  space.  The  height  of  the 
contour  surface  shows  the  value  of  the  progress  variable  in  the  flamelet  solutions,  while  the  color  contours 
show  the  value  of  temperature. 
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Once  Z  and  A  are  selected  as  independent  coordinates,  the  magnitude  of  their  associated 
fluxes  can  be  isolated  in  a  scalar  transport  equation  by  applying  flamelet-type  coordinate 
transformation  rules  such  as  [1], 

d(0  9Z  0(Q  dA  d(-)  de  d(-) 

dxj  dxj  dZ  dxj  dA  dxj  de 

This  particular  rule  is  valid  for  the  scalar  gradient  operator.  If  changes  along  the  third 
coordinate  e  are  assumed  to  be  small,  then  a  2-D  flamelet  equation  in  Z  and  A  space  can  be 
produced  using  this  transformation  [19].  After  the  neglect  of  unsteady  terms  and  the  use  of 
asymptotic  arguments,  the  transformed  equation  for  the  progress  variable  may  be  written, 

dO  v_  d^O 

M  [p^L)U|VA|  -  V  •  (pw A)]  -p^-—  =  pioc.  (10) 

The  first  set  of  terms  on  the  left  hand  side  of  this  equation  can  be  recognized  as  sharing 
the  form  of  the  transport  terms  in  the  premixed  flamelet  equations  (Eq.  (5)),  while  the 
second  term  on  the  left  hand  side  is  simply  the  diffusion  term  from  the  non-premixed 
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flamelet  equations  (Eq.  (2)).  The  right  hand  side  represents  the  chemical  source  term, 
demonstrating  again  that  the  local  combustion  regime  can  be  determined  by  examining 
which  mixing  process  acts  as  the  sink  for  the  chemical  source. 

A  model  for  determining  the  combustion  regime  at  a  given  temporal  and  spatial  location 
directly  follows  from  this  transformation.  The  terms  on  the  left  hand  side  of  Eq.  (10)  can 
be  grouped  according  to  the  regime  they  describe,  locally  calculated  in  the  flow  solver, 
and  then  compared  to  one  another  to  determine  each  regime’s  contribution  to  the  chemical 
source  term  budget.  More  explicitly,  the  terms  are  grouped  to  form  the  regime  specific 
quantities  0pre  and  0non , 


0 


pre 


dC 

=  [puSL,u  |  VA|  -  V  •  (pVV A)] , 


0, 


=  -p- 


Xz  ()2C 
2  dZ 2’ 


(11) 

(12) 


where  ©pre  describes  how  the  source  term  is  balanced  in  the  premixed  limit  and  0non 
describes  how  the  source  term  is  balanced  in  the  non- premixed  limit.  Note  that  outside 
of  the  premixed  flammability  limits,  the  laminar  burning  velocity  sl,u  effectively  becomes 
zero.  The  premixed  and  non-premixed  terms  can  be  very  simply  compared  to  determine 
their  relative  importance, 

q  _  ©pre 


0, 


(13) 


When  0  C  1,  the  regime  is  definitively  non-premixed.  Conversely,  the  regime  is  definitively 
premixed  when  0  1. 

The  0  variable  serves  to  characterize  the  local  combustion  regime  in  an  asymptotic 
sense,  but  it  does  not  describe  how  solutions  from  different  regimes  should  be  combined 
when  their  relative  contributions  to  the  chemical  source  term  are  similar  in  magnitude.  Such 
a  regime  combination  algorithm  is  the  second  task  required  of  a  multi-regime  model.  Here, 
a  simple  weighting  procedure  is  used  to  determine  the  local  contribution  of  each  regime. 
Flamelet  solutions  from  the  non-premixed  regime  will  be  labeled  as  4>k,non{Z,C ),  while 
flamelet  solutions  from  the  premixed  regime  will  be  labeled  <pk,pre(Z:  C).  These  solutions 
will  be  combined  at  any  temporal  and  spatial  location  according  to  the  rule 


0fc(£j  Z,  C)  —  4>k,non(Z ,  C)  ■  (1  —  £)  +  (f>k,pre{Z ,  C)  ■  £, 

where  the  weighting  coefficient  £  is 

Iv  dV 


£  = 


max(Jy  0pre  dV  fy  0 non  dV,  e) 


(14) 


(15) 


and  where  e  is  a  small,  positive  number  and  V  is  the  volume  of  a  domain  comprised  of  the 
neighboring  computational  cell  in  each  direction.  The  e  term  is  employed  to  ensure  that  the 
regime  is  always  well  defined,  while  the  integration  is  employed  to  ensure  that  transitions 
between  regimes  are  smooth.  The  limiting  cases  of  fully  premixed  and  fully  non-premixed 
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combustion  are  captured  by  this  weighting  procedure,  since  £  — >  1  when  Qpre  Qnon-,  and 

conversely  £  — >  0  when  Qnon  Qpre- 

It  is  important  to  note  that  the  quantities  A,  dC/d A,  and  d2C/dZ 2,  are  not  imme¬ 
diately  available  from  transport  equations.  These  quantities  can  either  be  solved  for  in 
flamelet  space  and  accessed  from  the  tabulated  flamelet  solutions  as  needed,  or  they  could 
be  calculated  by  considering  the  conditional  dependence  of  C  on  Z  and  A  within  a  region 
surrounding  the  spatial  location  of  interest.  In  the  flames  considered  here,  the  d2C /dZ2 
term  is  calculated  directly  from  the  non-premixed  flamelet  solutions,  and  stored  with  these 
solutions  so  that  it  can  be  accessed  using  the  Z  and  C  scalars.  Similarly,  A  and  dC/d  A 
could  be  calculated  from  non-premixed  flamelets  by  assigning  each  flamelet  an  appropriate 
value  of  the  A  variable  and  then  differentiating. 

When  quantities  such  as  d2C /dZ2  are  evaluated  using  non-premixed  flamelets,  a  mod¬ 
eling  assumption  is  being  invoked.  These  flamelet-based  evaluations  assume  that  non- 
premixed  transport  can  be  effectively  measured  by  treating  local  flame  structure  as  non- 
premixed.  If  this  non-premixed  treatment  is  valid,  then  it  is  of  course  valid  to  take  d2C /dZ2 
from  a  non-premixed  flamelet.  If  the  local  flame  structure  is  premixed,  however,  then 
d2C /dZ2  might  no  longer  be  accurately  described  by  a  non-premixed  flamelet.  Nonethe¬ 
less,  it  would  be  expected  that  in  a  premixed  combustion  regime  the  \z  term  that  is 
evaluated  from  the  flow  solver  would  be  small.  This  small  value  of  xz  would  tend  to  drive 
the  non-premixed  transport  term  Qnon  toward  zero,  and  minimize  the  influence  of  d2C /dZ2 . 
Additionally,  the  presence  of  significant  premixed  transport  would  drive  up  the  magnitude 
of  the  premixed  term  Qpre.  The  regime  indicator  would  accordingly  predict  premixed  com¬ 
bustion.  It  is  therefore  expected  that  d2C/dZ2  can  be  safely  evaluated  from  non-premixed 
flamelets  in  both  limiting  regimes.  The  laminar  flame  simulations  that  are  presented  below 
confirm  this  expectation. 

Although  the  A  variable  successfully  accounts  for  the  statistical  dependence  of  Z  and 
C,  it  is  recognized  to  be  a  challenging  quantity  from  both  a  physical  and  computational 
standpoint.  Furthermore,  its  definition  allows  for  some  uncertainty.  For  example,  A  is 
defined  in  the  context  of  non-premixed  flamelets,  but  Figs.  1-2  hint  that  a  possibility 
exists  of  finding  a  higher  value  of  C  in  premixed  flamelet  space  than  in  non-premixed 
flamelet  space.  If  this  difference  of  maximum  progress  variable  values  were  to  occur,  an 
arbitrary  means  of  extrapolating  a  value  of  A  to  a  higher  value  of  progress  variable  would 
have  to  be  introduced  in  premixed  regions  of  a  flow  field.  In  an  effort  to  minimize  these 
practical  challenges,  an  alternative  method  of  accounting  for  the  dependence  of  Z  and  C  is 
proposed. 

This  method  retains  the  use  of  a  flamelet  transformation  based  on  the  A  variable,  but 
once  the  Qpre  term  has  been  defined  transformation  rules  such  as  those  in  Eq.  (9)  are  again 
applied  to  Qpre-  These  transformation  rules  are  used  to  change  the  dependencies  on  A  to 
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dependencies  on  C  and  Z.  After  manipulation,  the  form  of  the  Opre  term  becomes 


©pre  — 


Pu  SL,h 


|VC|-^|VZ| 


v ' {pvvc)  -  <v ' {pvvz))  -pl^m 


(16) 


It  can  be  seen  that  this  method  subtracts  out  the  progress  variable’s  dependence  on  mixture 
fraction  in  an  explicit  fashion,  rather  than  through  a  specially  defined  variable.  Note  that 
the  ©non  term  remains  unchanged  relative  to  Eq.  (12)  in  this  modified  formulation.  In 
the  simulations  that  follow  the  regime  indicator  £  will  be  defined  according  to  Eq.  (21) 
using  Qpre  from  Eq.  (16)  and  Qnon  from  Eq.  (12).  Note  that  the  flame  simulations  shown 
below  in  sections  2.3.2  and  ??  indicate  that  a  regime  formulation  based  on  A  yields  results 
that  are  very  similar  to  the  results  of  the  modified  formulation.  Because  of  the  challenges 
associated  with  A,  the  results  of  the  more  tractable  modified  formulation  shown  in  Eq.  (16) 
are  presented  throughout  the  remaining  sections. 

Regardless  of  how  Qpre  is  calculated,  it  can  be  seen  that  both  sides  of  the  flamelet 
transformation  in  Eq.  (10)  approach  zero  in  non-reactive  areas  of  the  flow  field  [19].  The 
regime  indicator  £  then  ceases  to  provide  information.  This  does  not  cause  difficulties  in 
unburned  gas,  since  premixed  and  non-premixed  flamelets  have  identical  chemical  states  at  a 
given  mixture  fraction  and  zero  progress  variable.  At  fully  burned  conditions,  however,  this 
issue  needs  to  be  considered  because  the  production  of  pollutants  or  soot  may  be  sensitive  to 
the  composition  of  the  burned  gas.  Similarly,  regime  dependencies  need  to  be  considered  in 
the  preheat  regions  of  premixed  flames,  where  diffusion  leads  to  a  temperature  increase  but 
where  reaction  is  not  yet  significant  [19].  In  these  preheat  regions,  diffusive  and  convective 
premixed  transport  processes  cancel  with  one  another,  and  the  premixed  transport  term 
Qpre  would  diminish  to  zero  even  though  the  regime  is  definitively  premixed. 

These  issues  will  be  accounted  for  in  the  multi-regime  model  by  forcing  the  premixed 
term  Qpre  to  be  at  least  as  large  as  the  premixed  convective  flux  term, 
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This  bounding  operation  ensures  that  premixed  transport  will  be  appropriately  captured 
even  in  regions  where  premixed  diffusive  transport  and  premixed  convective  transport  op¬ 
pose  one  another.  Consequently,  premixed  transport  will  be  observed  even  in  preheat 
regions  of  premixed  flame  fronts,  where  chemical  source  terms  are  small.  Note  that  this 
bounding  operation  should  never  activate  in  reactive  gas,  because  in  the  presence  of  signif¬ 
icant  reaction  the  diffusive  term  is  expected  to  be  negative,  and  to  therefore  provide  a  net 
positive  addition  to  0pre. 

A  similar  kind  of  ambiguity  exists  in  the  non-premixed  flamelet  equations.  At  mixture 
fraction  values  away  from  stoichiometric  burning  conditions,  the  second  derivative  term 
in  Eq.  (12)  tends  toward  zero.  Consequently,  the  measured  non-premixed  transport  term 
©non  tends  toward  zero.  It  is  therefore  possible  to  find  flow  regions  that  are  non-premixed 
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iii  nature,  but  that  are  subject  to  only  very  weak  transport  in  mixture  fraction  space.  In 
an  effort  to  ensure  that  these  flow  regions  are  treated  correctly,  the  second  derivative  in  the 
non-premixed  term  Qnon  will  be  bounded  by  the  model  parameter  7, 

Qnon  =  ~P~2  '  mm  (  Q^2  >  7  )  ■  (18) 

Note  that  the  second  derivative  in  this  equation  is  taken  from  flamelet  space,  while  the 
dissipation  rate  is  taken  from  the  flow  solver.  Bounding  the  second  derivative  by  7  ensures 
that  the  regime  indicator  treats  marginally  reacting  regions  where  there  is  a  high  mixture 
fraction  dissipation  rate  and  almost  no  premixed  transport  as  non-premixed.  7  should  be 
much  smaller  than  the  value  of  the  second  derivative  at  the  most  reactive  mixture  fraction, 
so  that  it  does  not  influence  the  regime  model  in  reacting  flow  regions.  For  example,  in 
the  remaining  sections  the  7  parameter  is  set  equal  to  7  =  —  1,  which  is  less  than  1%  of 
the  peak  magnitude  of  the  second  derivatives  found  in  the  non-premixed  flamelets  shown 
in  section  2.2.1. 

Taken  together,  these  premixed  and  non-premixed  bounding  operations  serve  to  ensure 
that  the  local  regime  tends  toward  premixed  in  non-reactive  areas  where  xz  asymptotes  to 
zero,  and  towards  non-premixed  in  non-reactive  areas  where  \z  is  significant. 

2.2.3  Extending  The  Regime  Indicator  To  LES 

The  regime  indicator  can  be  extended  to  LES  by  considering  filtered  transport  terms.  For 
example,  after  LES  filtering  the  mixing  terms  that  result  from  the  coordinate  transformation 
can  be  written, 
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A  filtered  regime  indicator  may  then  be  formed  by  comparing  the  magnitudes  of  these 
mixing  terms, 


_ fy  Qpr  dV _ 

ma x(fv  0pr  dV  +  Jy  0np  dV,  e) 


(21) 


Once  the  regime  indicator  £  is  known,  the  premixed  flamelet  solutions  4>k,pr  and  the 
non-premixed  flamelet  solutions  (j>k,np  will  be  combined  for  use  in  a  flow  solver  according  to 
the  construct 


MC,  Z,  Z"2,C,  C”2)  =  (j>k,np(Z,  Z"2,  C)  ■  [1  -  £]  +  d>k, Pr(Z ,  c,  C"2)  ■  [ej  .  (22) 
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2.3  Validation  Of  The  Multi-Regime  Flamelet  Model 

In  this  section  a  laminar  flame  will  be  simulated  for  the  purpose  of  validating  the  multi¬ 
regime  flamelet  model.  The  flame  that  is  chosen  is  a  canonical  partially  premixed  triple 
flame.  This  flame  will  be  described  using  four  distinct  descriptions  of  chemistry:  1)  finite 
rate  chemistry;  2)  the  multi-regime  flamelet  model;  3)  a  purely  non-premixed  flamelet 
model;  and  4)  a  purely  premixed  flamelet  model. 

2.3.1  Triple  Flame  Description 

A  laminar  triple  flame  similar  to  the  flames  studied  by  Favier  and  Vervisch  [20]  and  Knud- 
sen  and  Pitsch  [19,  21]  is  selected  as  a  validation  case.  A  schematic  of  this  case  is  shown  in 
Fig.  4,  where  the  progress  variable  and  mixture  fraction  in  the  flame  is  plotted  at  initial¬ 
ization  and  then  again  at  a  time  when  the  flame  is  nearing  steady  state.  In  addition  to  the 
progress  variable,  Fig.  4  denotes  the  simulation’s  inlet  and  outlet,  as  well  as  the  mixture 
fraction  stratification  imposed  at  the  inlet.  All  of  the  triple  flame  simulations  are  initialized 
by  placing  two  counterrotating  vortices  in  the  middle  of  a  9  mm  x  6  mm  2-D  domain. 
These  vortices  are  introduced  in  order  to  force  the  flame  front  to  experience  a  variety  of 
transient  mixing  conditions.  Behind  the  vortices,  a  progress  variable  profile  in  the  form  of  a 
hyperbolic  tangent  along  the  downstream  direction  is  introduced.  The  mixture  fraction  in 
the  downstream  region  is  set  to  have  a  constant  stoichiometric  value.  Tabulated  premixed 
flamelet  solutions  are  then  accessed  using  the  Z  and  C  profiles  to  set  the  initial  values  of 
the  individual  chemical  species  and  temperature  in  the  domain.  The  inlet  of  the  simula¬ 
tion,  which  encompasses  the  entire  left  edge  of  the  domain,  is  set  to  have  a  constant  bulk 
value  of  downstream  velocity  (U=0.72  m/s)  and  zero  cross-stream  velocity.  The  species  and 
temperature  profiles  along  this  inlet  are  set  to  fully  unburnt  values,  but  a  mixture  fraction 
gradient  is  created  in  the  cross-stream  direction  along  the  inlet  so  that  the  effects  of  both 
premixed  and  non-premixed  combustion  can  be  considered.  The  mixture  fraction  at  the  top 
edge  of  the  inlet  is  set  as  .2=0.0,  while  the  mixture  fraction  at  the  bottom  edge  of  the  inlet 
is  set  as  Z= 0.5.  At  the  time  of  initialization,  this  cross-stream  gradient  of  mixture  fraction 
extends  downstream  to  the  vortices.  The  outlet  boundary  conditions,  which  are  active  on 
the  right  edge  of  the  domain,  are  solved  using  a  convective  equation.  The  boundaries  at 
the  top  and  bottom  of  the  domain  are  set  as  slip  walls.  The  simulation  is  run  until  the 
influence  of  the  vortices  has  disappeared  and  a  steady,  propagating  triple  flame  has  formed. 

2.3.2  Triple  Flame  Results 

Progress  Variable  Comparison 

Progress  variable  fields  from  the  finite  rate  triple  flame  solution  and  from  each  of  the 
flamelet  model  solutions  are  shown  in  Fig.  5.  The  plots  represent  information  from  three 
distinct  points  in  time,  with  time  increasing  from  left  to  right.  The  finite  rate  solution  is 
shown  in  the  top  row  of  the  figure,  while  the  flamelet  model  solutions  are  shown  in  the 
lower  rows. 
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t  1  Inlet 


Figure  4:  Schematic  of  the  triple  flame  simulation.  The  first  two  plots  show  the  progress  variable  C  and 
mixture  fraction  Z  at  initialization.  The  second  two  plots  show  these  quantities  once  the  propagating  triple 
flame  has  formed. 


The  flamelet  models  capture  the  general  flow  features  of  the  ‘true’  finite  rate  solution 
reasonably  well.  In  the  first  and  third  columns  of  Fig.  5,  for  example,  the  flamelet  models 
reproduce  all  of  the  global  structures  that  are  observed  in  the  finite  rate  case.  Although  these 
global  structures  are  accurately  described,  the  modeled  solutions  do  have  difficulty  capturing 
several  details  of  the  finite  rate  solution.  Especially  at  the  intermediate  time  shown  in  the 
figure’s  second  column,  the  flamelet  models  alter  the  evolution  of  flame-vortex  interactions 
relative  to  the  finite  rate  solution.  Additionally,  the  figure’s  third  column  shows  that  the 
propagation  speed  of  the  leading  triple  flame  edge  depends  on  the  modeling  approach.  These 
differences  in  the  progress  variable  field  hint  at  the  more  significant  differences  that  appear 
when  details  of  the  flame  structure  are  considered. 

Regime  Predictions 

The  regime  indicator  £  from  Eq.  (21)  is  plotted  in  the  multi-regime  model  simulations 
in  Fig.  6.  In  this  plot,  black  coloring  (£  =  1)  indicates  that  the  local  combustion  regime 
is  premixed,  while  white  coloring  (£  =  0)  indicates  the  local  regime  is  non-premixed.  The 
indicator  predicts  that  the  leading  edge  of  the  flame  front  is  definitively  premixed,  while  the 
trailing  region  of  the  flame  tends  to  be  largely  non-premixed.  These  results  are  consistent 
with  the  sharp  progress  variable  gradients  that  occur  at  the  flame  front,  and  with  the  weak 
mixture  fraction  diffusion  that  occurs  across  the  middle  of  the  burnt  gas  behind  the  flame 
front. 
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Time  (s|  0.000900 
Multi-Regime 


Time  |s|  0.000901 
Premixed 


Time  |s|  0.000900 
Non-Preinixed 


Figure  5:  Progress  variable  fields  in  the  triple  flame,  computed  from  the  finite  rate  approach  and  from  the 
multi-regime,  the  purely  premixed,  and  the  purely  non-premixed  flamelet  models. 


Figure  6:  Regime  indicator  £  from  Eq.  (21)  in  the  multi-regime  triple  flame  simulation.  £  =  1  (black) 
denotes  that  premixed  flamelet  solutions  are  accessed,  while  £  =  0  (white)  denotes  that  non-premixed 
flamelet  solutions  are  accessed. 


Edge  Speed  Comparison 

The  influence  of  the  combustion  regime  can  be  analyzed  in  more  detail  by  considering  the 
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Figure  7:  Position  and  speed  of  the  leading  flame  edge  in  the  triple  flame  from:  finite  rate  chemistry  (black 
•  );  multi-regime  model  (red  ■);  premixed  model  only  (green  ♦);  non-premixed  model  only  (blue  A).  The 
position  of  the  leading  flame  edge  is  normalized  by  the  length  of  the  computational  domain,  Lq. 


position  and  the  speed  of  the  leading  flame  edge.  These  quantities  are  plotted  in  Fig.  7  for 
each  of  the  triple  flame  cases.  The  leading  flame  edge  location  is  defined  as  the  horizontal 
coordinate  in  the  computational  domain  where  the  progress  variable  field  first  reaches  a 
value  of  (7=0.08.  This  leading  edge  is  largely  premixed  in  nature,  as  described  by  both 
the  multi-regime  model  and  the  flame  index.  Consequently,  the  multi-regime  and  premixed 
models  describe  the  evolution  of  the  leading  edge  more  accurately  than  the  non-premixed 
model.  The  non-premixed  model’s  tendency  is  to  under-predict  the  true  flame  speed. 

Figure  8  shows  a  sampling  of  the  conditional  source  term  values  that  are  responsible 
for  the  flame  propagation  speeds  in  Fig.  7.  Note  that  the  premixed  and  multi-regime 
source  terms  are  taken  from  stored  premixed  flamelets,  but  are  dependent  on  the  level  set 
as  described  in  [14].  The  non-premixed  model’s  under-prediction  of  the  flame  edge  speed 
is  due  to  the  under-predicted  source  terms  shown  in  Fig.  8.  The  largest  source  terms  in 
the  flame  are  found  around  a  mixture  fraction  of  Z=0.075,  and  the  non-premixed  sources 
under-predict  the  finite  rate  sources  at  this  mixture  fraction.  The  linear  region  of  the  non- 
premixed  source  term  profile,  located  between  C=0.0  and  C=0.16  in  the  top  two  plots 
of  Fig.  8,  describes  the  interpolated  non-premixed  flamelet  space  that  sits  between  the 
unburned  mixing  line  and  the  burning  flamelet  with  the  lowest  maximum  temperature. 

At  a  much  richer  mixture  fraction  of  Z=0.14,  Fig.  8  demonstrates  that  the  non-premixed 
model  alters  its  relative  behavior  and  over-predicts  the  finite  rate  source  term.  Because 
sources  at  rich  Z  values  are  relatively  small,  however,  this  over-prediction  does  not  sig¬ 
nificantly  influence  flame  propagation  in  the  non-premixed  model.  Just  as  at  the  leaner 
mixture  fractions,  source  terms  along  the  Z=0.14  surface  are  most  accurately  described  by 
the  multi-regime  and  premixed  models. 

A  final  noteworthy  feature  of  the  edge  speed  plot  in  Fig.  7  is  the  slight  difference  that 
exists  between  the  premixed  and  the  multi-regime  predictions.  This  difference  is  not  due 
to  the  progress  variable  source  term  that  is  used  for  modeling,  since  Fig.  8  demonstrates 
that  the  premixed  and  multi-regime  sources  are  virtually  identical.  The  difference  is  instead 
a  product  of  how  premixed  flamelet  solutions  are  extrapolated  beyond  their  flammability 
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Progress  Variable,  C  Progress  Variable,  C 


Figure  8:  Conditional  progress  variable  source  terms  in  the  triple  flame,  from:  finite  rate  chemistry  simulation 
(black  •);  multi-regime  model  (red  ■);  premixed  model  only  (green  ♦);  non-premixed  model  only  (blue  A). 


limits.  Figure  9  shows  typical  density  fields  in  the  finite  rate,  multi-regime,  and  premixed 
simulations.  In  these  density  fields  fuel-rich,  unburned  gas  is  seen  to  flow  around  the  rich 
(lower)  edge  of  the  triple  flame.  As  shown  on  the  right  in  Fig.  4,  the  mixture  fraction 
associated  with  this  rich,  unburned  gas  is  Z= 0.5.  This  Z  value  lies  outside  the  premixed 
flammability  limits,  and  the  density  at  (Z= 0.5,  C=0.0)  in  the  premixed  model  is  therefore 
calculated  by  linearly  interpolating  between  the  unburned  density  at  the  rich  flammability 
limit  and  the  density  at  a  mixture  fraction  of  unity.  Because  density  actually  varies  non- 
linearly  between  Z= 0  and  Z= 1,  this  interpolation  introduces  an  error.  The  error  is  evident 
when  the  premixed  density  plot  in  Fig.  9  is  compared  to  the  finite  rate  density  plot.  The 
premixed  model’s  over-prediction  of  density  affects  the  deflection  of  velocity  streamlines, 
and  consequently  slows  down  the  leading  edge  speed.  Density  in  the  premixed  model  could 
of  course  be  corrected  using  a  different  extrapolation  procedure,  but  some  ambiguity  would 
nevertheless  remain. 

Flame  Structure  Comparison 

Plots  showing  the  evolution  of  the  CO  species  in  the  triple  flame  simulations  are  provided 
in  Fig.  10.  Just  as  in  the  progress  variable  plots,  some  model  dependencies  appear  during 
the  transient  flame  development  phase  shown  in  the  middle  column.  For  example,  a  strip  of 
high  CO  concentration  is  seen  along  the  simulation  centerline  at  the  right  edge  of  the  domain 
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Figure  9:  Density  fields  in  the  triple  flame  simulations,  computed  from  the  finite  rate  approach  (left)  and 
from  the  multi-regime  (center)  and  purely  premixed  (right)  models. 


in  the  finite  rate  simulation.  This  high  CO  concentration  is  reproduced  by  the  multi-regime 
model,  but  is  absent  in  both  the  purely  premixed  and  the  purely  non-premixed  results. 
Further  dependencies  appear  in  the  third  column  of  Fig.  10,  especially  on  the  fuel-rich 
(lower)  side  of  the  triple  flame.  At  this  later  time  the  multi-regime  and  premixed  models 
again  reasonably  reproduce  the  finite  rate  case,  with  the  exception  of  an  over-prediction  of 
CO  in  the  rich  trailing  flame  edge.  Conversely,  the  non-premixed  model  under-predicts  CO 
throughout  the  rich  side  of  the  flame. 

Conditional  species  and  temperature  data  from  two  different  mixture  fraction  values  are 
plotted  in  Figs.  11-12.  Unlike  the  progress  variable  source  terms,  species  and  temperature 
data  are  only  a  mild  function  of  the  regime  at  Z=0.075.  Fig.  11  shows  that  the  non- 
premixed  model  somewhat  under-predicts  the  conversion  of  CO  to  CO2,  and  also  slightly 
under-predicts  the  OH  concentration  at  low  progress  variable.  Although  these  differences 
are  small,  the  premixed  and  multi-regime  models  agree  more  closely  with  the  finite  rate 
case  than  does  the  non-premixed  model. 

Conditional  data  from  the  richer  mixture  fraction  Z=0.14  is  plotted  in  Fig.  12.  Here 
the  comparisons  are  more  interesting,  and  the  species’  dependence  on  the  local  combustion 
regime  is  more  significant.  At  first  glance,  the  finite  rate  CO  profile  in  Fig.  12  appears  to  be 
most  accurately  predicted  by  the  purely  non-premixed  model.  For  example,  at  intermediate 
C  values  the  overlapping  multi-regime  and  premixed  results  deviate  from  the  finite  rate  case 
much  more  than  the  non-premixed  model  results. 

An  incorrectly  predicted  regime  does  not  explain  the  entirety  of  Fig.  12,  however.  Three 
observations  suggest  that  the  errors  in  the  multi-regime  CO  predictions  are  due  to  more 
than  just  the  regime  indicator.  First,  the  non-premixed  model  is  seen  to  over-predict  CO 
oxidation  in  the  burned  gas.  This  over-prediction  is  indicated  by  the  sharp  downturn  in 
the  non-premixed  CO  profile  in  Fig.  12  around  (7=0.22.  Second,  the  lower  half  of  Fig.  12 
demonstrates  that  the  non-premixed  model  does  not  predict  temperature  or  OH  any  more 
accurately  than  the  premixed  model.  Third,  the  upper  right  plot  of  Fig.  1  shows  that,  at 
the  mixture  fraction  Z=0.14  considered  in  Fig.  12,  the  lowest  CO  value  that  can  be  found  in 
a  non-premixed  burning  flamelet  is  Y co  =  0.06.  Consequently,  between  progress  variable 
values  of  (7=0.0  (where  Y co  =  0.0)  and  (7=0.15  (where  Y co  =  0.06),  the  non-premixed  CO 
profile  in  Fig.  12  is  the  result  of  a  simple  linear  interpolation.  The  non-premixed  model’s 
agreement  with  finite  rate  chemistry  in  this  region  is  consequently  due  less  to  physical 
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Figure  10:  Vco  fields  in  the  triple  flame,  computed  from  the  finite  rate  approach  and  from  the  multi-regime, 
premixed,  and  non-premixed  flamelet  models. 


modeling  accuracy  than  to  numerical  good  fortune. 

Therefore,  while  the  CO  profile  in  Fig.  12  does  suggest  that  the  regime  is  not  fully 
premixed,  the  finite  rate  data  could  neither  be  reproduced  by  a  fully  non-premixed  model. 
It  appears,  then,  that  chemistry  on  the  fuel-rich  side  of  the  flame  deviates  from  1-D  flamelet 
manifolds.  Heat  diffusion  across  mixture  fraction  surfaces  is  one  possible  source  of  these 
deviations.  If  a  fuel-rich  leading  edge  of  the  triple  flame  were  premixed,  for  example,  heat 
diffusion  from  an  adjacent  and  higher  temperature  premixed  reaction  zone  might  cause  the 
local  temperature  to  rise  above  the  relevant  1-D  premixed  flamelet  temperature.  This  kind 
of  behavior  is  indeed  observed  in  the  temperature  plot  in  Fig.  12.  Accounting  for  higher 
order  flame  structure  effects  [22]  such  as  the  diffusion  processes  associated  with  stratified 
flame  propagation  might  improve  the  finite  rate  and  multi-regime  model  agreement. 

When  comparing  flame  structure,  it  should  be  emphasized  that  the  non-premixed  regime 
must  be  included  in  the  model  formulation  to  obtain  the  best  agreement  with  the  finite  rate 
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Figure  11:  Flame  structure  profiles  in  the  triple  flame  conditioned  on  Z=0.075  from:  finite  rate  chemistry 
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simulation.  This  was  first  demonstrated  in  the  more  accurate  flame  edge  speed  predictions 
of  the  multi-regime  model  relative  to  the  purely  premixed  model.  It  was  also  demonstrated 
by  the  multi-regime  model’s  ability  to  most  accurately  represent  many  of  the  finite  rate 
flame  structures  that  appear  during  flame  transients,  as  shown  in  the  middle  column  of 
Fig.  5. 

NO  Comparison 

Figure  13  shows  the  evolution  of  the  NO  species  in  the  triple  flame  simulations.  Again, 
the  finite  rate  simulation  is  shown  in  the  first  row,  followed  by  the  various  flamelet  models. 
The  third  column  of  Fig  13  shows  that  while  the  premixed  model  comes  close,  none  of 
the  flamelet  models  predict  as  much  NO  production  as  the  finite  rate  simulation.  This 
is  consistent  with  the  presence  of  strain  and  stratification  within  the  triple  flame,  and 
with  the  observation  that  the  finite  rate  flame  temperatures  were  slightly  higher  than  the 
premixed  flamelet  model  temperatures.  To  further  investigate  the  regime  dependency  of 
NO  production,  conditional  plots  of  both  NO  and  the  production  rate  of  NO  are  shown  in 
Fig.  14. 

The  data  on  the  left  side  of  Fig.  14  is  taken  at  an  early  time  in  the  simulation,  while 
the  data  on  the  right  side  is  taken  at  the  later  time  associated  with  the  third  column  of 
Fig.  13.  The  Z=0.075  surface  that  is  used  for  conditioning  is  very  near  to  the  mixture 
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(black  •);  multi-regime  model  (red  ■);  premixed  model  only  (green  ♦);  non-premixed  model  only  (blue  A). 


fraction  value  at  which  NO  production  reaches  its  maximum.  Although  many  of  the  multi¬ 
regime  flamelet  predictions  up  to  this  point  have  closely  corresponded  to  the  premixed 
flamelet  predictions,  Fig.  14  indicates  that  the  opposite  situation  occurs  with  NO.  The 
multi-regime  predictions  closely  match  those  of  the  non-premixed  flamelet  model  because 
NO  forms  over  relatively  long  time  scales,  and  the  regime  indicators  in  Fig.  6  show  that 
regions  far  downstream  of  the  flame  front  are  predicted  to  be  non-premixed.  Interestingly, 
these  regime  observations  have  little  practical  consequence  for  the  triple  flame,  since  the 
NO  production  rate  does  not  significantly  vary  between  regimes.  The  similarity  of  these  NO 
production  rates  is  consistent  with  the  underlying  flamelet  solutions:  the  non-premixed  NO 
production  rate  shown  in  Fig.  1  and  the  premixed  NO  production  rate  shown  in  Fig.  2  both 
reached  maximums  of  approximately  0.005  kg  /  m3  s.  Note,  however,  that  NO  production 
is  strongly  a  function  of  scalar  dissipation  rate  (see  Fig.  1),  and  in  general  NO  production 
may  differ  significantly  across  combustion  regimes. 

2.3.3  Summary 

Within  the  context  of  the  flamelet  modeling  approach,  an  increasing  need  exists  for  models 
and  validation  data  that  are  relevant  to  multi-regime  combustion  processes.  Sections  2.2 
and  2.3  of  this  report  attempt  to  address  a  part  of  that  need  by  considering  1)  whether 
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Figure  13:  NO  fields  in  the  triple  flame,  computed  from  the  finite  rate  approach  and  from  the  purely 
non-premixed,  the  purely  premixed,  and  the  multi-regime  flamelet  models. 


a  mixed  regime  flamelet  approach  can  capture  combustion  and  pollutant  formation  in  a 
detailed  chemistry,  multi-dimensional  environment,  and  2)  whether  the  errors  associated 
with  the  flamelet  model  are  attributable  to  incorrect  regime  descriptions,  or  to  problems 
with  the  underlying  single-regime  flamelet  assumptions. 

To  address  these  issues,  simulations  of  a  primarily  premixed  triple  flame  were  performed 
using  finite  rate  chemistry  and  a  variety  of  flamelet  combustion  models.  It  was  demonstrated 
that  a  multi-regime  flamelet  approach  could  capture  the  correct  asymptotic  combustion 
regime  in  the  flame.  Further  details  of  this  study  can  be  found  in  references  [19]  and  [4], 
which  were  supported  by  this  project. 

Although  the  multi-regime  approach  correctly  identified  combustion  regimes,  several  of 
the  corresponding  flamelet  model  predictions  nonetheless  deviated  from  the  finite  rate  data. 
These  deviations  are  most  likely  due  to  strain  and  stratification.  In  the  flamelet  modeling 
context,  these  deviations  can  only  be  accounted  for  by  expanding  the  parameters  used  to 
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Figure  14:  NO  and  NO  production  source  term  profiles  in  the  triple  flame  conditioned  on  Z=0.075,  and  at 
two  different  times:  finite  rate  chemistry  (black  •);  multi-regime  model  (red  ■);  premixed  model  only  (green 
♦);  non-premixed  model  only  (blue  ▲). 


describe  chemistry.  In  spite  of  these  modeling  limitations,  the  multi-regime  approach  was 
shown  to  accurately  distinguish  between  the  premixed  and  non-premixed  asymptotes.  When 
the  triple  flame  and  the  counterflow  flame  are  considered  as  a  single  modeling  challenge, 
the  advantage  of  the  multi-regime  model  is  that  it  can  be  applied  to  both  cases  without 
the  need  for  case-related  mixing  assumptions.  The  ability  to  apply  a  single  model  that 
dynamically  accounts  for  regime  information  led  to  significantly  improved  accuracy  relative 
to  the  single  regime  approaches. 

2.4  Modeling  Sub-Filter  Scalar  Dissipation  In  Turbulent  Non-Premixed 
Flames 

2.4.1  Motivation 

Scalar  dissipation  rates  and  sub-filter  scalar  variances  are  critical  modeling  parameters  in 
large  eddy  simulations  (LES)  of  reacting  flows.  A  variety  of  models  for  these  quantities 
have  been  proposed,  many  having  been  designed  to  operate  within  a  Reynolds  averaged 
framework.  Dissipation  rate  and  variance  models  can  be  categorized  according  to  whether 
they  are  based  on  an  algebraic  equation  or  a  transport  equation.  Algebraic  models  for 
the  variance  [23,  24,  25,  26,  27]  and  the  dissipation  rate  [23,  28,  29]  directly  describe  the 
modeled  parameter  using  information  about  the  local  scalar  field  and  filter  width,  and  have 
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the  advantage  of  being  conceptually  simple  and  computationally  inexpensive.  Transport 
equation  models  for  the  variance  [27,  30,  31]  and  for  the  dissipation  rate  [32,  33,  34]  offer 
different  advantages:  they  do  not  forcibly  assume  that  production  and  dissipation  processes 
are  in  equilibrium,  they  produce  fields  that  are  less  noisy  than  algebraic  models,  and  they 
incorporate  a  wider  range  of  physics.  Nonetheless,  transport  equation  approaches  are  sen¬ 
sitive  to  the  descriptions  used  to  model  unclosed  sub-filter  terms,  and  the  performance  of 
these  equations  is  not  always  superior  to  the  performance  of  algebraic  models  [27].  For  ex¬ 
ample,  the  closure  of  these  transport  equations  may  be  strongly  affected  by  the  details  and 
accuracy  of  the  numerical  method  that  is  used.  Consequently,  many  uncertainties  continue 
to  surround  the  modeling  of  scalar  variance  and  dissipation. 

Two  observations  regarding  the  sub-filter  dissipation  motivate  this  section.  First,  im¬ 
provements  in  variance  and  dissipation  rate  LES  models  are  expected  to  significantly  im¬ 
prove  the  quality  of  reactive  flow  predictions.  Second,  significant  ambiguity  exists  regarding 
how  these  improvements  should  be  pursued.  The  goal  of  this  work  is  to  further  address 
these  modeling  needs  by  making  use  of  newly  available  reacting  direct  numerical  simulation 
data,  and  by  proposing  a  new  LES  closure  model. 

2.4.2  Modeling  Background 

The  scalar  dissipation  rate  xz  of  a  scalar  Z  whose  diffusivity  is  T>z  can  be  written 

Xz  =  2Vz\S7Z\2.  (23) 

When  an  LES  filter  is  applied  to  this  quantity,  the  filtered  dissipation  rate  is  found, 

Xz  =  2  VZ  |VZ|2,  (24) 

where  the  operator  denotes  spatial  filtering.  The  sub-filter  scalar  variance  is  typically 
calculated  using  density  weighted  filtering  that  will  be  denoted  by  the  tilde  operator  ( ~ ) . 
This  sub-filter  variance  is 

=  Z2-  (Z)2.  (25) 

Algebraic  Models 

Algebraic  LES  models  for  \z  can  be  formulated  using  a  variety  of  parameters  [29].  The 
simplest  and  most  widely  used  formulation  separates  the  dissipation  rate  into  resolved  and 
sub-filter  components  [28], 

Xz  =  ‘2{Vz  +  Vt)\VZ\2,  (26) 

where  the  filtered  molecular  diffusivity  T>z  is  associated  with  the  resolved  dissipation,  and 
the  turbulent  diffusivity  T>t  is  associated  with  the  sub-filter  dissipation.  This  model  in¬ 
troduces  no  unknown  coefficients,  as  the  turbulent  diffusivity  is  already  available  from  the 
solution  of  the  underlying  scalar’s  transport  equation. 

Algebraic  LES  models  for  the  sub-filter  scalar  variance  can  also  take  many  forms  [24, 
25,  26].  The  most  widely  used  form  relates  the  variance  to  gradients  of  the  resolved  scalar 
field, 

Zffs  =  Cvar  A2  jVZ\2.  (27) 
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The  A  parameter  in  Eq.  (27)  is  the  LES  filter  width.  This  model  introduces  a  unique 
coefficient  that  can  also  be  dynamically  calculated  using  a  variety  of  methods  [10,  26]. 


Transport  Equation  Models 

An  exact  transport  equation  can  be  written  for  the  scalar  dissipation  rate.  Similarly, 
the  scalar  variance  can  be  described  exactly  using  the  solution  of  either  a  variance  equation 
or  transport  equations  for  both  Z  and  Z 2 .  In  the  context  of  LES,  several  of  the  terms  that 
appear  in  the  equations  for  these  quantities  require  closure  models.  Here,  exact  equations 
for  the  dissipation  rate  and  the  square  of  a  scalar  (Z2)  are  introduced  first. 

In  reacting  flows,  the  transport  equation  for  xz  includes  terms  that  depend  on  gradients 
of  diffusivity  and  density.  When  this  equation  is  solved  in  LES,  these  terms  all  require 
closure.  To  limit  the  appearance  of  unclosed  terms,  a  transport  equation  for  |VZ|2,  rather 
than  for  xz,  is  considered.  Letting  D/Dt  denote  the  material  derivative,  this  equation  is 
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(28) 


(a) 

(b) 

(c) 

(d) 

(e) 


The  (a)  and  ( b )  terms  in  Eq.  (28)  describe  the  production  and  dissipation,  respectively,  of 
the  scalar  dissipation  rate.  The  (c)  through  (e)  terms  describe  the  effects  of  changes  in 
density  and  diffusivity,  and  reduce  to  zero  when  p  and  T>z  are  constant.  Once  solved,  the 
|VZ|2  quantity  from  Eq.  (28)  can  be  multiplied  with  T>z  to  determine  xz- 

Applying  an  LES  density  weighted  filter  to  Eq.  (28)  results  in  a  transport  equation  for 

the  term  \S7Z\2.  It  will  be  assumed  that  correlations  between  mixture  fraction  gradients 
and  the  scalar  diffusivity  are  small,  so  that  the  transported  quantity  can  be  multiplied  with 
the  filtered  diffusivity  to  compute  the  dissipation  rate, 


2VZ  |VZ|2  «  2Vz\VZ\2. 


(29) 


A  priori  tests  have  suggested  that  the  density  weighting  which  differentiates  the  filtering  of 
jVZj2  and  |  VZ|2  can  be  neglected  on  the  grounds  that  it  has  a  relatively  small  influence  [35]. 


31 


After  filtering,  the  equation  for  the  quantity  |VZ|2  becomes, 


wt^z i2)  =^(^  +  D‘^(^2)) 
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n  ^  f  d2Z  \2 
P  z  \dxidxj 

2  dp  dZ  7  d  f  dZ  \  \ 

p  dxi  dxi  \  dxj  v  Z  dxj  )  J 

,  MPVZ)  dZ  (  921T 
dxi  dxi  \dxjdxj  J 

l2dz  dz  fd2(pvz)\ 

dxi  dxj  \  dxidxj  J 

where  a  standard  turbulent  diffusivity  model  has  been  invoked  to  describe  sub-filter  scalar 
flux.  The  (a)  through  (e)  terms  in  Eq.  (30)  are  unclosed.  Sanders  et  al.  [34]  provides  a 
review  of  the  relevant  closure  modeling  for  Reynolds  averaged  approaches.  In  LES,  different 
closures  are  needed.  This  issue  will  be  further  addressed  in  section  2.4.5. 

The  sub-filter  scalar  variance  can  be  solved  for  in  a  number  of  different  ways.  Based  on 
Kaul  et  aV s  observations  [27],  a  transport  equation  for  Z 2  will  be  used  here.  After  LES 
filtering,  this  transport  equation  is  written 

w  i^2) = £]  + <3i) 

Equation  (31)  is  noteworthy  in  that  it  does  not  add  to  the  closure  modeling  problem. 
Rather,  if  the  scalar  dissipation  rate  \z  has  been  closed,  then  Eq.  (31)  is  closed.  Once 
solved,  Eq.  (31)  can  be  used  to  calculate  the  sub-filter  variance  Z”js  using  Eq.  (25). 

2.4.3  Auto-Igniting  Jet 

DNS  Description 

The  DNS  case  that  will  be  used  to  assess  the  dissipation  and  variance  models  is  an 
auto-igniting  slot  jet.  This  ethylene  fueled  flame  was  originally  the  subject  of  a  DNS  study 
by  Yoo  et  al.  [36].  A  schematic  of  the  flame  is  shown  in  Fig.  15,  where  velocity  contours 
are  plotted  on  the  left  and  temperature  contours  are  plotted  on  the  right.  The  central  fuel 
jet  has  a  bulk  velocity  of  220  m/s,  a  temperature  of  550  K,  and  a  composition  of  82%  N2 
and  18%  C2H4  by  volume.  This  central  jet  stream  is  denoted  by  the  mixture  fraction 
composition  Z  =  1.  The  co-flow  that  surrounds  the  central  jet  consists  of  100%  air  at  a 
temperature  of  1550  K.  The  air  co-flow  enters  the  domain  at  a  bulk  velocity  of  20  m/s,  and 


(a) 

(b) 

(c) 

(d) 

(e) 
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this  air  stream  is  denoted  by  a  mixture  fraction  of  Z  =  0.  No  nozzle  separates  the  main 
jet  and  co-flow;  rather,  the  velocity,  temperature,  and  composition  profiles  are  prescribed 
to  smoothly  transition  between  Z  =  0  and  Z  =  1  according  to  the  functions  listed  in 
reference  [36].  The  total  width  of  the  central  jet  is  H  =  0.002  m,  and  the  jet  Reynolds 
number  associated  with  the  DNS  conditions  is  10,000. 


Chemistry  in  the  flame  is 
described  using  a  22  species 
ethlyene  mechanism  [36]  that 
was  reduced  from  a  larger  hy¬ 
drocarbon  mechanism.  This 
same  22  species  mechanism  is 
used  in  all  modeled  LES  calcu¬ 
lations. 

LES  Description 

In  all  LES  simulations  of 
the  jet  DNS,  combustion  is 
described  using  an  unsteady 
flamelet  approach  similar  to 
that  of  Ihme  and  See  [37]  and 
Pitsch  [38] .  The  model  is  imple¬ 
mented  by  solving  the  unsteady 
flamelet  equations  [1,  12].  Unsteady  flamelet  solutions  are  generated  for  a  variety  of  refer¬ 
ence  scalar  dissipation  rates.  They  are  then  tabulated  as  a  function  of  Z.  Xz.ref ,  and  the 
progress  variable  C. 

Once  generated,  the  unsteady  solutions  are  convoluted  with  presumed  PDFs  for  appli¬ 
cation  in  the  LES.  A  beta  PDF  is  presumed  to  describe  Z .  while  a  delta  PDF  is  presumed 
to  describe  Xz,ref  ■  Additionally,  it  is  assumed  that  a  single  unsteady  flamelet  solution  is 
representative  of  the  conditions  in  an  LES  mesh  cell.  Consequently,  only  the  mean  value 
of  the  progress  variable  is  needed  for  modeling  [10].  These  assumptions  lead  to  a  tabulated 
LES  chemistry  database  of  the  form 


Figure  15:  The  ethylene  flame  DNS.  Left:  axial  velocity  contours; 
Right:  temperature  contours. 


fa  =  MZ1Z"2,XZ,re{,C), 


(32) 


where  </>£  is  any  reacting  quantity  of  interest.  The  reference  dissipation  rate  can  be  deter¬ 
mined  from  the  local  unconditional  LES  dissipation  rate  Xz  by  filtering  Eq.  (3), 


XZ, ref  =  Xz  ' 


f(zre{) 

m  ’ 


(33) 


where  f(Z)  is  determined  by  convoluting  the  function  f(Z)  with  a  beta  PDF. 

Filtered  C  and  Z  scalars  are  transported  in  the  LES,  and  the  scalar  dissipation  rate  and 
sub-filter  variance  are  calculated  using  models  from  section  2.4.2.  The  tabulated  chemistry 
database  from  Eq.  (32)  is  accessed  during  the  LES  to  determine  chemical  information  such 
as  the  density  and  the  progress  variable  source  term. 
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Two  LES  runs  are  performed  us¬ 
ing  different  mesh  resolutions.  The 
first  mesh  consists  of  1  million  (1M) 
cells,  and  corresponds  to  an  LES  Li¬ 
ter  width  to  Kolmogorov  scale  ratio 
of  approximately  A/77  =  8-  The  sec¬ 
ond  LES  mesh  consists  of  23  million 
(23M)  cells,  and  the  corresponding 
length  scale  ratio  is  A/77  =  4-5.  The 
DNS  data  is  filtered  for  comparison 
with  LES  using  coarse  and  fine  fil¬ 
ter  widths  that  correspond  to  the  1 
million  cell  and  23  million  cell  LES 
meshes,  respectively. 

2.4.4  Algebraic  Model  Perfor¬ 
mance 

Figure  16  presents  fields  from  the 
1  million  cell  and  23  million  cell 
LES  runs  alongside  DNS  fields  that 
are  filtered  using  the  A  from  the  1 
million  cell  LES.  The  filtered  mix¬ 
ture  fraction  is  plotted  along  with 
the  algebraically  modeled  dissipa¬ 
tion  rate  from  Eq.  (26)  and  the  al¬ 
gebraically  modeled  sub- filter  jjcalar 
variance  from  Eq.  (27).  The  Z  field 
in  the  higher  resolution  LES  is  seen 
to  capture  significantly  more  flow 
structure  than  the  lower  resolution 
LES.  Consequently,  the  burden  on 
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Figure  16:  LES  algebraic  model  results.  Left  column:  DNS 
fields  filtered  using  the  1  million  cell  LES  filter  width;  Middle 
column:  1  million  cell  LES;  Right  column:  23  million  cell  LES. 
Upper  row:  Z  contours;  Middle  row:  \z  fr°m  Eq.  (26);  Lower 
row:  Z"js  from  Eq.  (27). 


the  sub-filter  models  is  eased  in  the 
higher  resolution  data  set.  The  1 
million  cell  LES  is  therefore  viewed 
as  the  primary  sub-filter  modeling 
challenge. 

Figure  17  compares  the  alge¬ 
braically  modeled  scalar  dissipation  rate  from  Eq.  (26)  with  the  filtered  dissipation  rate 
from  the  DNS.  All  quantities  are  time  averaged,  as  represented  by  the  (•)  operator.  Fig¬ 
ure  17(a)  shows  Xz  as  a  function  of  the  cross-stream  coordinate  at  various  downstream 
planes,  while  Fig.  17(b)  shows  the  mixture  fraction  conditioned  Xz  iR  those  same  planes. 
The  algebraic  model’s  representation  of  Xz  is  shown  to  be  poor,  especially  in  the  lower 
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(a)  physical  space  (b)  mixture  fraction  space 

Figure  17:  Algebraic  \z  m°del  results  in  the  reacting  jet,  plotted  in  both  physical  (a)  and  conditional  (b) 
space.  DNS  (•••);  1  million  cell  LES  ( - );  23  million  cell  LES  ( - ). 


resolution  case.  At  both  LES  resolutions,  the  model  under-predicts  \z  from  the  DNS.  In 
the  lower  resolution  case  in  particular,  the  under-predictions  are  order-of-magnitude  dif¬ 
ferences.  Because  ignition  is  sensitive  to  local  dissipation,  these  errors  lead  to  significant 
under-predictions  of  the  flame  lift-off  height.  Although  the  Xz  predictions  are  improved 
when  the  LES  resolution  is  increased,  differences  between  the  high  and  low  resolution  LES 
indicate  that  the  algebraic  model  is  undesirably  sensitive  to  the  filter  size. 

Resolved  and  sub-filter  contributions  to  the  LES  algebraic  Xz  model  are  separated  and 
compared  with  the  DNS  in  Fig.  18.  Both  LES  filter  widths  are  considered.  The  resolved 
component  of  dissipation  is  calculated  as  (2D^|VZ|2)  in  both  the  filtered  DNS  and  LES, 
while  the  sub-filter  component  is  calculated  as  {2T>z\ VZ|2  —  2T>z\ VZ|2)  and  (2Dt|VZ|2)  in 
the  DNS  and  LES,  respectively.  Figure  18(a)  plots  the  dissipation  rate  components  that 
are  associated  with  the  coarse  filter  width.  At  this  filter  resolution,  the  resolved  dissipation 
from  the  filtered  DNS  (solid  black  circles)  and  sub-filter  dissipation  from  the  filtered  DNS 
(open  blue  circles)  are  comparable  in  magnitude,  sub-filter  dissipation  modeling  therefore 
plays  an  important  role  in  the  1  million  cell  LES,  even  though  the  Reynolds  number  of  the 
case  is  not  representative  of  realistic  combustor  conditions.  The  resolved  DNS  dissipation 
rate  is  reasonably  well  described  by  the  resolved  component  of  the  LES  model  (black  solid 
line),  at  least  in  the  rich  part  of  the  flame  where  turbulence  is  most  intense.  The  sub¬ 
filter  LES  dissipation  (blue  dashed  line),  however,  universally  underestimates  the  sub-filter 
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Figure  18:  Budget  of  the  algebraic  \z  m°del  in  the  reacting  jet.  Resolved  dissipation  from  the  filtered  DNS 

(•••);  resolved  dissipation  from  the  algebraic  LES  model  ( - );  sub- filter  dissipation  from  the  filtered 

DNS  (o  o  o);  sub-filter  dissipation  from  the  algebraic  LES  model  ( - ).  The  filter  resolutions  correspond 

to  the  1  million  cell  LES  (a)  and  the  23  million  cell  LES  (b). 


dissipation  seen  in  the  DNS.  This  sub- filter  modeling  error  is  responsible  for  the  bulk  of  the 
error  in  the  total  modeled  dissipation  rate. 

Figure  18(b)  demonstrates  that  the  sub-filter  component  of  the  algebraic  model  performs 
poorly  even  when  the  filter  width  is  decreased.  In  spite  of  this  poor  performance,  the 
total  dissipation  rate  is  predicted  more  accurately  (see  Fig.  17).  The  explanation  for  the 
improvement  is  that  the  dissipation  rate  budget  shifts  from  the  sub-filter  scales  to  the 
resolved  scales.  The  sub-filter  model  therefore  reduces  in  importance.  Dissipation  rate 
predictions  in  the  23  million  cell  LES  can  be  said  to  improve  relative  to  the  1  million  cell 
LES  predictions  in  spite  of,  and  not  because  of,  sub-filter  model  performance. 

2.4.5  An  Adapted  Dynamic  LES  Closure  For  The  Dissipation  Equation 

One  method  of  addressing  the  observed  sub-filter  dissipation  errors  would  be  to  adjust  the 
form  of  the  algebraic  model  [25,  29].  Here,  however,  the  transport  equations  that  were 
shown  in  section  2.4.2  are  considered  as  a  means  of  improving  the  modeling.  This  section 
addresses  the  closure  of  the  LES  transport  equation  for  |VZ|2,  shown  in  Eq.  (30). 

Production  And  Dissipation  Closure  Assumptions 

Five  unclosed  terms  labeled  (a)  through  (e)  appear  in  the  dissipation  rate  equation, 
Eq.  (30).  The  two  most  important  of  these  are  the  production  and  dissipation  terms.  They 
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are  split  into  resolved  and  sub- filter  components  for  use  in  LES, 
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Following  Reynolds  averaged  modeling  approaches  such  as  those  reviewed  by  Sanders  [34, 
39],  the  sub- filter  production  and  dissipation  terms  will  be  written  for  LES  as 


(  '3  \V2  _ 

Vsfs  =  Cprd-p-  -(|VZ|2-|VZ|2), 

nsfs  =  Cdis-p-SZ-  (l  VZ|2  -  |VZ|2)2 . 

Zsfs 


(36) 

(37) 


Cprd  and  Cdis  are  model  coefficients,  and  Z”2s  is  the  sub- filter  scalar  variance  that  is  modeled 
either  algebraically  or  with  a  transport  equation,  vt  is  the  turbulent  viscosity  that,  in  the 
following  sections,  will  be  calculated  using  a  dynamic  Smagorinsky  model  with  Lagrangian 
averaging  [40].  The  riA/ A  term  in  the  production  model  can  be  viewed  as  the  sub-filter 
momentum  dissipation,  ea-  In  the  homogeneous  DNS  calculations  presented  next,  is 
locally  calculated  using  the  sub-filter  kinetic  energy  ksfs, A  as  n'A  =  ((2/3)/cs/Sja)1//2-  In 
the  LES  calculations  of  the  reacting  jet  that  are  shown  later,  nA/ A  is  calculated  following 
Deardorff  [41,  42], 

It  will  be  assumed  that  the  remaining  (c)  through  (e)  terms  in  Eq.  (30)  are  of  less 
importance  on  sub-filter  scales,  as  justified  in  [42], 


Modeled  LES  Scalar  Dissipation  Equation 

These  dissipation  closure  assumptions  lead  to  the  following  form  of  the  |VZ|2  equation, 
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Once  closed,  |VZ|2  from  Eq.  (38)  can  be  multiplied  with  T>z  to  determine  xz-  This 
Xz  value  can  then  be  used  as  the  source  term  in  the  sub- filter  variance  equation,  Eq.  (31). 
The  only  remaining  barrier  to  implementation  of  the  transport  equation  models  for  the 
dissipation  and  variance  is  the  specification  of  Cpr(i  and  Cdis- 


Scalar  Mixing  Turbulence  Database 

Closure  models  for  Cprd  and  Cdis  are  especially  difficult  to  develop  because  most  of 
the  dissipation  rate’s  energetic  content  exists  on  the  smallest  scales  of  a  flow.  This  small 
scale  content  inhibits  the  accuracy  of  traditional  dynamic  LES  approaches  for  describing 
model  coefficients,  which  tend  to  extrapolate  information  from  large  resolved  scales  to  small 
unresolved  scales.  When  the  large  scales  contain  negligible  amounts  of  information  relative 
to  the  small  scales,  the  extrapolation  procedure  ceases  to  aid  in  modeling.  Alternatives  to 
the  traditional  dynamic  algorithm  are  therefore  needed. 

The  evaluation  of  Cprd  and  Cdis  will  be  investigated  using  a  second  database  of  di¬ 
rect  numerical  simulations  that  describe  conserved  scalar  mixing  in  the  presence  of  forced, 
constant  density,  homogeneous  isotropic  turbulence.  This  second  set  of  DNS  data  was  orig¬ 
inally  developed  in  order  to  study  the  probability  distribution  functions  that  result  from 
multi-scalar  mixing  [43,  44].  These  DNS  cases  are  run  on  a  5123  computational  mesh  and 
are  forced  using  the  scheme  of  Rosales  and  Meneveau  [45]  to  a  Reynolds  number  of  approxi¬ 
mately  Rey=100.  Scalar  fields  are  then  initialized  in  the  domain  according  to  the  procedure 
of  Eswaran  et  al.  [46].  Further  details  regarding  the  scalar  mixing  DNS  can  be  found  in  [44]. 

Data  is  extracted  from  the 
homogeneous  DNS  for  analysis 
at  one  eddy  turnover  time  af¬ 
ter  the  scalar  field  is  initial¬ 
ized.  Typical  contour  planes 
that  show  the  scalar  field  at 
this  point  in  time  appear  in 
Fig.  19.  This  scalar  data  is 
filtered  for  analysis  in  an  LES 
context  using  clipped  Gaussian 
filter  kernels  [14].  The  filtering 
is  performed  using  several  filter 
widths,  and  model  results  are 
Figure  19:  Typical  planes  from  two  of  the  the  5123  scalar  mixing  ho-  parameterized  by  these  widths, 
mogeneous  turbulence  DNS  runs  [43],  showing  contours  of  the  scalar. 


Model  Coefficient  PDFs 

Figure  20  shows  the  PDFs 

that  describe  the  distribution  of  the  model  coefficient  Cprd  after  one  eddy  turnover  time  in 
the  scalar  mixing  DNS.  This  coefficient  is  evaluated  at  each  location  in  the  homogeneous 
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prd 


Figure  20:  PDFs  of  the  model  coefficient  Cprd  from  Eq.  (38),  evaluated  from  the  homogeneous  scalar  mixing 
DNS.  The  PDFs  are  parameterized  by  the  LES  filter  width,  A. 


DNS  domain  using  the  formula 


Cprd  — 


-(z/jA)1/2 


az  az  _  o-aa  dz  oz 

r  dxi  dxi  dxj  r  dxi  dxi  dxd 


p-K)3/2.(|VZ|2-|VZ|- 


(39) 


PDFs  of  the  coefficient  are  constructed  for  a  given  filter  size  by  binning  results  from  the 
entire  domain. 

Figure  20  shows  that  the  Cprd  PDF  is  uni-modal.  Additionally,  the  value  of  the  model 
coefficient  is  reasonably  filter  independent:  PDFs  associated  with  different  filter  widths 
have  approximately  the  same  mean,  and  relax  towards  a  single  value  as  the  filter  width 
increases.  Cprd  is  therefore  reasonably  represented  by  a  single,  filter-independent  constant. 
In  agreement  with  several  Reynolds  averaged  modeling  approaches  [34],  Fig.  20  indicates 
that  this  coefficient  should  have  a  value  of  unity. 

PDFs  of  the  model  coefficient  Cfns  are  shown  in  Fig.  21,  where  C^is  is  calculated 
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The  PDFs  in  Fig.  21  are  again  unimodal,  but  the  value  of  Cdis  is  filter  dependent.  Conse¬ 
quently,  this  coefficient  cannot  be  described  as  a  single  value.  Filter  dependencies  such  as 
these  are  often  dealt  with  by  employing  dynamic  models,  but  the  dynamic  procedure  is  not 
applicable  to  quantities  like  |VZ|2  whose  energy  content  exists  on  small  scales. 


Adapted  Dynamic  LES  Closure 

Here  an  alternative  to  the  traditional  dynamic  approach  is  proposed.  This  alternative 
is  rooted  in  the  recognition  that  many  turbulence  modeling  coefficients  respond  similarly 
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Figure  21:  PDFs  of  Cdia  from  Eq.  (38),  evaluated  from  the  homogeneous  scalar  mixing  DNS.  The  PDFs  are 
parameterized  by  the  LES  filter  width  A. 


to  changes  in  the  local  intensity  of  turbulence.  For  example,  if  the  passage  of  a  turbulent 
eddy  increases  the  local  turbulent  scalar  flux,  the  scalar  variance  might  be  expected  to 
increase  in  the  same  way.  Similarly,  a  flow  event  that  causes  the  local  scalar  variance  to 
increase  might  also  increase  the  local  dissipation  rate.  Viewed  from  the  reverse  perspective, 
sub-filter  model  coefficients  should  tend  toward  zero  in  flow  regions  that  laminarize  or  that 
are  adjacent  to  walls. 

The  idea  of  the  adapted  dynamic  closure  is  to  take  information  from  quantities  that 
are  amenable  to  dynamic  calculations,  and  to  apply  that  information  to  other  quantities 
that  are  not  similarly  amenable.  This  projection  of  information  is  useful  because  small  scale 
quantities  are,  like  large  scale  quantities,  sensitive  to  local  turbulence.  Information  regarding 
local  turbulence  can  be  used  to  increase  the  accuracy  of  a  model  coefficient  regardless  of 
where  the  modeled  quantity’s  spectrum  peaks. 

This  projection  of  information  is  possible  because  an  appropriate  value  of  the  model 
coefficient  must  exist.  Indeed,  the  problem  with  applying  a  standard  dynamic  procedure 
to  small  scale  quantities  is  not  that  the  coefficients  are  indeterminate.  Rather,  the  prob¬ 
lem  is  that  coefficient  values  cannot  be  calculated  due  to  the  absence  of  useful  large  scale 
information.  The  adapted  dynamic  closure  simply  circumvents  this  calculation  problem  by 
borrowing  information  from  other  readily  obtained  dynamic  solutions. 

The  adapted  dynamic  closure  can  be  used  in  the  scalar  mixing  DNS  to  compute  the 
coefficient  C(ils .  Dynamic  information  will  be  borrowed  from  the  algebraic  scalar  variance 
model  in  Eq.  (27),  and  C(i%s  will  be  rewritten  as, 


Cdis  =  CdiStd- 

V  ^  var  / 


(41) 


Cdis,d  in  this  expression  is  now  the  unknown  model  coefficient,  while  Cvar  is  the  dynami¬ 
cally  computed  variance  coefficient  from  Eq.  (27).  C®ar  is  a  baseline  value  of  the  variance 
coefficient  that  can  be  determined  by  using  a  Taylor  series  expansion  to  derive  the  algebraic 
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Figure  22:  PDFs  of  the  model  coefficient  Cdis,d  from  Eq.  (41),  evaluated  from  the  homogeneous  scalar  mixing 
DNS.  The  PDFs  are  parameterized  by  the  LES  filter  width,  A. 


variance  model  [26].  Proceeding  through  this  derivation  leads  to  C®ar  =  1/12.  This  number 
can  be  viewed  as  the  variance  coefficient’s  value  in  a  particular  turbulence  regime.  The  de¬ 
tails  of  this  regime,  and  of  its  relationship  to  homogeneous  turbulence,  are  irrelevant.  These 
details  will  be  accounted  for  in  the  specification  of  the  constant  Cdis,d ■  The  intensity  of  Cvar 
relative  to  this  referenced  regime  is  the  only  information  of  interest,  and  this  information  is 
captured  by  the  ratio  in  Eq.  (41). 

The  adapted  dynamic  closure  is  tested  by  using  the  scalar  mixing  DNS  to  calculate  PDFs 
of  the  coefficient  Cdis,d ■  Coefficients  are  determined  by  first  using  the  DNS  to  determine 
the  exact  value  of  Cvar  at  all  locations  in  the  DNS  domain.  CdiS,d  is  then  determined  from 
the  expression 
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where  (7°or  =  1/12.  The  sub-filter  variance  Z"2S  in  this  expression  is  computed  directly 
from  the  DNS,  and  the  coefficient  Cvar  is  computed  so  that  the  algebraic  sub-filter  variance 
model  is  exact.  Equation  (42)  differs  from  the  expression  used  to  calculate  C!rils  only  in  that 
it  is  multiplied  by  the  inverse  of  the  ratio  ( Cvar  /  C°ar). 

PDFs  describing  the  distribution  of  Cdis,d  in  the  scalar  mixing  DNS  are  plotted  in  Fig.  22. 
These  distributions  are  largely  independent  of  the  LES  filter  width,  and  are  unimodal.  The 
information  that  is  provided  by  the  dynamic  variance  calculation  therefore  does  account  for 
the  filter  dependencies  of  the  Cdis  coefficient.  The  removal  of  this  information  reduces  the 
extent  of  the  physics  that  a  model  coefficient  must  attempt  to  describe.  When  the  adapted 
dynamic  procedure  is  incorporated  into  the  modeling  framework,  the  remaining  coefficient 
Cdis,d  can  be  accurately  treated  as  a  single,  constant  value. 
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The  mean  values  of  the  CdiS,d  distributions  are  approximately  CdiS,d  =  2.3,  but  the  peaks 
of  the  distributions  occur  around  CdiS,d  =  0.75.  Many  RANS  models  use  a  dissipation  model 
coefficient  of  approximately  1.0  [34]  that  sits  in  between  these  peak  and  mean  values.  To 
complete  the  closures  for  the  LES  |VZ|2  equation  the  value  Cdis,d  =  1-0  is  used,  and  the 
coefficient  CdiS  in  Eq.  (38)  is  modeled  as 

Cdis  =  Cdis,d  •  ^  =  1.0  •  =  12  Cvar.  (43) 

The  LES  transport  equation  for  |VZ|2  may  then  be  written  in  final  form  as 
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where  Cprd= 1.0  and  Cvar  is  the  dynamically  computed  sub-filter  scalar  variance  coefficient 
from  Eq.  (27).  The  variance  Z”2S  that  appears  in  the  last  term  in  Eq.  (44)  is  determined 
from  the  solution  of  the  conserved  Z  equation  and  the  Z2  equation  (Eq.  (31)). 


2.4.6  Transport  Equation  Model  Performance 

The  transport  equation  models  for  \z  an(l  ^sfs  are  now  evaluated  by  returning  to  the  LES 
simulations  of  the  reactive  jet  case.  Additional  computations  of  both  the  1  million  cell  LES 
and  the  23  million  cell  LES  are  performed,  and  the  |VZ|2  and  Z 2  variables  are  solved  for 
using  Eqs.  (44)  and  (31),  respectively.  These  two  variables  are  used  to  locally  compute  the 
dissipation  rate  Xz  and  the  sub- filter  variance  Z”jg. 

Time  averaged  transported  dissipation  rates  from  both  LES  simulations  are  compared 
to  the  time  averaged  DNS  dissipation  rate  and  to  the  original  algebraic  model  results  in 
Fig.  23.  The  transported  dissipation  rate  model  improves  upon  the  algebraic  model  in  two 
ways.  First,  mesh  dependencies  are  reduced.  The  algebraic  model  results  were  shown  to 
be  a  strong  function  of  the  LES  filter  width.  Conversely,  the  transport  equation  results  in 
Fig.  23  are  only  weakly  dependent  on  the  LES  filter.  This  change  represents  a  significant 
improvement  in  model  robustness.  Second,  the  transported  dissipation  rate  in  the  1  million 
cell  LES  is  in  significantly  better  agreement  with  the  DNS  than  the  1  million  cell  LES 
algebraic  model. 
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(a)  1  million  cell  LES  filter  width  (b)  23  million  cell  LES  filter  width 
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(c)  1  million  cell  LES  filter  width  (d)  23  million  cell  LES  filter  width 


Figure  23:  Dissipation  rate,  m  both  physical  and  mixture-fraction  conditioned  space  in  the  reacting  jet. 

DNS  (•  •  •);  algebraic  \z  model  ( - );  transported  \z  model  ( - )•  (a)  and  (c):  1  million  cell  LES 

filter  width;  (b)  and  (d):  23  million  cell  LES  filter  width. 
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The  improved  agreement  is  attributed  to  two  factors:  1)  to  improved  prediction  of  the 
sub-filter  variance  production  and  variance  dissipation  processes  that  were  inadequately 
represented  by  the  algebraic  approach,  and  2)  to  the  relaxation  of  the  assumption  that  the 
dissipation  rate  is  always  in  equilibrium  with  the  local  rate  of  variance  production. 

2.4.7  Summary 

This  work  has  used  an  auto- igniting  jet  DNS  to  analyze  the  performance  of  LES  algebraic 
and  transport  equation  models  for  both  the  scalar  dissipation  rate  and  the  sub-filter  scalar 
variance.  It  was  shown  that  the  algebraic  models  under-predicted  both  the  dissipation 
and  the  sub-filter  variance  in  the  jet.  The  under-predictions  were  almost  entirely  due  to 
the  sub-filter  components  of  the  algebraic  expressions.  These  model  errors  motivated  an 
examination  of  the  closure  of  a  transported  scalar  dissipation  rate  equation.  An  additional 
DNS  of  scalar  mixing  in  homogeneous  turbulence  was  introduced  to  aid  this  examination. 
The  scalar  mixing  DNS  demonstrated  that  sub- filter  production  of  \z  could  be  described 
using  a  single  model  coefficient,  while  the  coefficient  associated  with  sub- filter  dissipation 
of  Xz  was  dependent  on  the  LES  filter.  To  account  for  this  filter  dependence,  an  adapted 
dynamic  procedure  was  introduced  and  applied  to  the  dissipation  model  coefficient.  This 
procedure  used  information  from  a  dynamic  algebraic  variance  model  to  reduce  the  amount 
of  physics  that  the  dissipation  model  must  describe.  The  resulting  transport  equation 
closure  for  Xz  was  shown  to  perform  well  for  a  variety  of  LES  filter  widths.  The  transport 
equation  models  were  then  applied  in  the  reacting  jet  LES.  The  transported  Xz  model 
performed  significantly  better  than  the  algebraic  model  and  improved  the  agreement  with 
the  DNS  results. 

2.5  Modeling  Sub-Filter  Strain  Effects  In  Turbulent  Premixed  Flames 
2.5.1  Motivation 

Premixed  combustion  is  difficult  to  describe  using  Large  Eddy  Simulation  (LES)  because 
premixed  flame  structures  are  typically  thinner  than  computationally  tractable  LES  filter 
widths.  The  interaction  of  mixing  physics  and  chemical  kinetics  is  therefore  under-resolved. 
Separation  between  the  resolved  LES  scales  and  the  flame  structure  scales  creates  mod¬ 
eling  challenges  when  mixing  is  strong  enough  to  perturb  flame  structures.  For  example, 
premixed  flame  structures  respond  to  the  influence  of  differential  diffusion  when  non-unity 
Lewis  number  fuels  are  used,  and  to  the  influence  of  turbulence  when  the  flame  is  in  the 
thin  reaction  zones  or  broken  reaction  zones  regime  of  the  premixed  regime  diagram  [1], 
Quantities  of  interest  such  as  burning  velocities  and  pollutant  formation  rates  can  be  highly 
sensitive  to  the  details  of  these  flame  structure  perturbations.  Modeling  these  perturbations 
is  therefore  critically  important  in  premixed  LES. 

Describing  flame  structure  perturbations  continues  to  be  a  critical  modeling  challenge, 
and  the  question  of  how  to  best  describe  these  perturbations  in  LES  remains  open.  In  this 
section  a  two  coordinate  ffamelet  model  is  applied  in  an  a  posteriori  LES  of  a  highly  strained 
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jet  flame  to  address  this  challenge.  A  tractable  implementation  of  a  flamelet  model  that 
accounts  for  strain  will  be  proposed  and  then  validated  using  a  DNS  database. 

2.5.2  Case  Description:  Premixed  Jet  DNS 

The  flame  case  that  is  modeled  is  a  pre¬ 
mixed  methane  and  air  slot-jet  direct  nu¬ 
merical  simulation  (DNS)  from  Sandia  Na¬ 
tional  Laboratories  [47,  48].  A  schematic 
of  this  case  is  presented  in  Fig.  24.  The 
DNS  was  run  at  atmospheric  pressure  us¬ 
ing  a  17  species  chemical  mechanism,  with 
4  of  these  species  assumed  to  be  in  steady 
state  [47].  The  velocity  of  the  central  jet 
stream  is  100  m/s  and  the  velocity  of  the 
coflow  stream  is  25  m/s.  The  mixture  frac¬ 
tion  of  the  two  streams  is  identical  and  cor¬ 
responds  to  an  equivalence  ratio  of  cj)  =  0.7. 
The  composition  in  the  central  jet  is  set 
as  unburned  CH4  and  air  at  a  temperature 
of  Tu  =  800  K,  while  the  composition  of 
the  coflow  is  set  as  the  combustion  prod¬ 
ucts  that  result  from  solving  a  premixed 
unstrained  flame  with  unburned  conditions 
equal  to  those  in  the  central  jet.  The  lami¬ 
nar  burning  velocity  and  flame  width  associ¬ 
ated  with  these  conditions  are  sl  =  1.8  m/s 
and  Ip  =  0.3  mm,  respectively.  The  width 
of  the  central  jet  is  H  =  1.8  mm  and  the  corresponding  jet  Reynolds  number  is  Re,  =  2100. 

The  magnitude  of  the  turbulent  fluctuations  in  the  central  jet  is  v! /sl  =  10  and  the 
integral  turbulent  length  scale  is  It /If  =  4.  The  corresponding  Karlovitz  (Ka  =  1‘p/rj2) 
and  Damkohler  (Da  =  SLlt/u'lp)  numbers  are  Ka  =  225  and  Da  =  0.4,  implying  that 
this  flame  sits  just  inside  the  broken  reaction  zones  regime  of  the  turbulent  premixed  flame 
regime  diagram  [1,  49].  Turbulence  is  therefore  expected  to  strongly  influence  premixed 
flame  structures.  The  DNS  was  run  on  a  195  million  (1200  x  600  x  270)  cell  mesh. 

2.5.3  Strained  Premixed  Flamelet  Model 

Strained  Flamelet  Formulation 

This  section  introduces  the  strained  premixed  flamelet  LES  model.  The  model  builds 
upon  laminar  strained  flamelet  models  [22,  50,  51,  52]  and  turbulent  unstrained  flamelet 
models  [19,  49,  53]  by  considering  strain  in  the  context  of  LES.  Chemistry  in  the  model 
will  be  described  using  solutions  of  the  counterflow  premixed  flamelet  equations  that  were 


Figure  24:  Premixed  CH4  jet  DNS  from  references  [47, 
48].  Left:  instantaneous  C  =  0.065  isosurface.  Right: 
density  weighted  and  time  averaged  C  field,  ranging 
from  0.00  (blue)  to  0.19  (red),  with  the  C  =  0.065 
isosurface  denoted  as  a  black  line. 
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developed  by  Dixon-Lewis  et  al.  [54],  These  equations  describe  premixed  flames  that  burn 
in  a  back-to-back  fashion  while  subject  to  strain.  Because  the  full  form  of  the  equations  is 
available  in  reference  [54],  they  are  only  briefly  reviewed  here.  The  equations  are  expressed 
as  a  function  of  a  1-D  similarity  coordinate  y  defined  as 


V  = 
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(45) 


The  y  variable  is  a  physical  space  coordinate  that  is  aligned  perpendicular  to  the  counterflow 
flame  fronts  and  parallel  to  the  direction  of  the  incoming  unburned  gas  flow.  The  density  is 
p,  the  viscosity  is  v,  and  the  subscript  0  denotes  a  reference  point  in  the  unburned  gas.  A 
coordinate  x  is  aligned  perpendicular  to  the  y  coordinate,  and  the  strain  rate  in  the  flame 
system  is  defined  as  the  derivative  of  the  velocity  in  the  x  direction  (u),  evaluated  in  the 
unburned  gas,  a  =  duo/dx. 

The  similarity  coordinate  can  be  used  to  write  species  equations  for  the  asymptotic 
counterflow  premixed  flames  as  a  function  of  one  variable.  For  example,  the  transport 
equation  for  the  species  4>i  is 

-/  d(Pi  -1  d  rhi 
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where  ViiV  is  the  standard  diffusion  velocity  of  4h  written  in  terms  of  the  y  direction.  The 
function  /  appears  in  expressions  for  the  physical  space  velocities  u  and  v  and  is  defined  so 
that  continuity  is  satisfied, 


pv  =  -p0f  ,  pu  =  xp0(df/dy).  (47) 

Inserting  these  velocities  into  the  momentum  transport  equations  yields  an  expression  that 
can  be  solved  to  determine  /.  Further  details  about  the  strained  equations  can  be  found  in 
reference  [54], 

In  the  present  study,  the  flamelet  equations  will  be  used  to  describe  a  symmetric  coun¬ 
terflow  flame  configuration,  and  will  therefore  be  solved  only  between  one  inlet  stream  and 
the  symmetry  plane.  Boundary  conditions  in  the  inlet  stream  are  set  using  information 
from  the  central  jet  in  the  DNS:  unburned  CH4  and  air  at  an  equivalence  ratio  of  0  =  0.7 
and  a  temperature  of  Tu  =  800K.  Neumann  conditions  are  applied  at  the  symmetry  plane. 

The  FlameMaster  program  [16]  is  used  to  solve  the  stretched  flamelet  equations  (Eq.  (46), 
e.g.)  with  the  17  species  CH4  mechanism  that  was  employed  in  the  DNS.  Differential  dif¬ 
fusion  is  expected  to  be  important  within  the  inner  reaction  zone,  and  constant  non-unity 
Lewis  numbers  are  therefore  used  for  all  chemical  species.  Flamelet  solutions  are  gener¬ 
ated  for  a  variety  of  values  of  the  strain  rate  parameter  a,  and  the  standard  s-curve  [1] 
that  results  from  these  calculations  is  shown  in  of  Fig.  25(a).  This  figure  shows  that  the 
maximum  flamelet  temperature  will  initially  decrease  as  the  imposed  strain  is  increased.  A 
turning  point  is  eventually  reached,  however.  Increasing  the  strain  beyond  the  critical  value 
associated  with  this  turning  point  has  the  effect  of  quenching  the  flamelet.  Decreasing  the 
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Figure  25:  Premixed  strained  flamelets.  (a):  maximum  flamelet  temperature  vs.  imposed  strain  rate  (the 

so-called  s-curve).  (b):  progress  variable  source  terms  in  unstrained  (•  •  •)  and  strained  ( - )  flamelet 

solutions.  The  rightmost  ends  of  the  strained  flamelet  solutions  in  (b)  correspond  to  the  symmetry  boundary 
condition,  and  the  strained  solutions  represent  the  entire  range  of  strainrates  shown  in  (a). 


strain  imposed  on  a  flamelet  near  the  turning  point,  conversely,  may  push  the  solution  onto 
the  so-called  middle  branch  of  the  s-curve.  This  branch  describes  steady  flamelet  solutions 
which  are  unstable  in  the  sense  that  they  will  relax  toward  the  upper  or  lower  branches  when 
perturbed.  Although  unstable,  these  middle  branch  solutions  are  important  in  the  premixed 
slot-jet  DNS  where  turbulence  induces  unsteadiness  in  the  premixed  flame  structure. 

A  progress  variable  will  be  defined  in  this  study  as  C  =  Yh2o  +  Yh2  +  Yco2  +  Yxo, 
where  the  Tj’s  denote  species  mass  fractions.  Figure  25(b)  shows  the  chemical  source  term 
associated  with  C  in  several  of  the  strained  flamelet  solutions.  Movement  down  the  s-curve 
decreases  the  magnitude  of  these  source  terms.  This  demonstrated  sensitivity  to  strain  is  a 
leading  order  effect  in  the  flamelet  model,  and  it  will  be  shown  below  that  this  effect  must 
be  captured  if  the  slot-jet  DNS  is  to  be  predicted  correctly. 

Tabulation  Coordinate  Selection 

The  stretched  flamelet  solutions  are  parameterized  and  tabulated  for  use  in  LES.  Changes 
along  a  single  flamelet  solution  are  parameterized  using  the  progress  variable  C.  A  second 
coordinate  is  needed  to  parameterize  the  influence  of  the  imposed  strain  rate  a.  The  dissi¬ 
pation  rate  of  the  progress  variable,  x.C  =  2I?c'|VCj2,  the  mass  fraction  of  a  minor  species, 
and  the  strain  rate  itself  could  all  be  used  as  a  second  coordinate.  It  is  important  to  note 
that  spatial  variations  of  each  of  these  coordinates  are  characterized  by  length  scales  on 
the  order  of  the  premixed  flame  structure,  and  that  this  presents  modeling  challenges.  Any 
coordinate  that  parameterizes  strain  will  be  subject  to  this  challenge,  however,  because 
coordinates  that  only  exist  on  larger  scales  will  not  be  sensitive  to  turbulent  strain  effects. 
The  particular  choice  of  a  second  coordinate  that  is  used  in  this  study  will  be  motivated  by 
referring  to  profiles  from  the  strained  flamelet  solutions  plotted  in  Fig.  26. 
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Figure  26:  (a):  progress  variable  dissipation  rate  in  an  unstrained  flamelet  (•••),  and  in  strained  flamelets 

above  ( - )  and  below  (— - )  the  turning  point  of  the  s-curve.  (b):  hydrogen  mass  fractions  in  the 

strained  flamelets. 


Figure  26(a)  shows  xc  as  a  function  of  C  in  several  flamelets.  The  limiting  solution  of 
zero  strain  (a  — >  0)  is  shown  as  solid  circles,  while  solutions  associated  with  the  upper  branch 
of  the  s-curve  in  Fig.  25(a)  are  shown  as  solid  lines.  On  the  upper  branch,  the  dissipation 
rate  tends  to  increase  as  the  imposed  strain  rate  increases.  Dissipation  peaks  near  the 
turning  point  of  the  s-curve  and  then  begins  to  fall  as  the  strain  rate  is  decreased.  Flamelet 
solutions  from  the  middle  branch  of  the  s-curve  are  shown  as  dashed  lines  in  Fig.  26(a). 
Many  of  these  middle  branch  solutions  are  characterized  by  dissipation  rates  that  are  similar 
to  the  rates  in  the  upper  branch  solutions.  The  \C  parameter  therefore  fails  to  uniquely 
parameterize  strain:  both  upper  and  middle  branch  solutions  can  be  found  at  a  single 
location  in  ( C ,  xc)  space.  This  lack  of  uniqueness  is  analogous  to  the  lack  of  uniqueness 
associated  with  parameterizing  non-premixed  flamelet  solutions  using  the  mixture  fraction 
and  its  dissipation  rate  [10]. 

Unlike  the  dissipation  rate,  the  hydrogen  radical  mass  fraction  ( Yu )  does  uniquely  pa¬ 
rameterize  the  influence  of  strain.  Figure  26(b)  shows  Yjj  as  a  function  of  C  in  flamelets 
spanning  the  entire  s-curve.  Movement  down  the  s-curve  results  in  a  monotonic  decrease 
of  hydrogen  over  virtually  all  of  progress  variable  space.  Additionally,  the  hydrogen  radi¬ 
cal’s  concentration  varies  significantly  as  a  function  of  the  imposed  strain.  This  variation 
is  helpful  in  that  small  numerical  errors  made  in  the  calculation  of  Yjj  are  not  expected 
to  influence  the  flamelet  solution  that  is  accessed.  The  hydrogen  mass  fraction  is  therefore 
selected  as  a  second  parameterizing  coordinate,  and  the  stretched  flamelet  solutions  are 
tabulated  as 

<j>i  =  <j>i(C,YH).  (48) 

The  use  of  hydrogen  captures  the  same  kind  of  information  as  the  coordinates  suggested 
by  other  researchers  in  fully  resolved  contexts.  For  example,  strain  effects  have  also  been 
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described  using  elemental  mass  fractions  [51],  the  CO  radical  [55],  the  OH  radical  [22], 
and  the  temperature  [52],  The  advantages  of  hydrogen  include  its  sensitivity  to  strain, 
its  straightforward  transport  equation,  and  its  small  Lewis  number  which  allows  for  some 
implicit  capture  of  differential  diffusion  effects.  Additionally,  the  use  of  a  chemical  species 
in  place  of  a  strain  rate  is  expected  to  reduce  the  sensitivity  of  the  model  to  the  details  of 
how  strain  is  induced  in  the  flamelet  equations  (back-to-back  flamelets  versus  a  different 
configuration) . 

Unstrained  Flamelet  Model  For  Comparison 

The  strained  flamelet  model  will  be  compared  with  an  unstrained  model  to  better  under¬ 
stand  the  importance  of  strain  effects  in  the  jet  flame  DNS.  In  place  of  the  strained  flamelet 
equations  (Eq.  (46),  e.g.),  this  unstrained  model  solves  1-D  unstrained  premixed  flamelets. 
This  model  is  representative  of  a  typical  premixed  flamelet  LES  implementation  [53,  19]. 
The  equation  that  governs  the  species  mass  fraction  in  the  unstrained  flamelets  is 

d 

PuSL'udx 

The  diffusion  velocity  in  the  x  direction,  Vj.x,  accounts  for  typical  molar  molecular  diffusion 
effects.  The  chemical  source  term  is  ci^,  the  unstrained  laminar  burning  velocity  is  sl,u , 
and  the  density  in  the  unburned  gas  is  pu.  Just  as  in  the  strained  flamelet  model,  constant 
non- unity  Lewis  numbers  are  used  for  all  species  and  the  equations  are  solved  with  the  17 
species  CH4  mechanism  employed  in  the  DNS. 

Because  no  strain  parameter  appears  in  the  unstrained  model,  only  one  flamelet  can 
be  used  to  describe  the  conditions  in  the  DNS  ( 4>  =  0.7,  Tu  =  800K).  The  resulting 
flamelet  solution  is  parameterized  using  the  progress  variable  C,  and  flamelet  quantities  in 
the  unstrained  model  are  accessed  as 
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Sub-filter  Modeling 

LES  filtering  will  be  accounted  for  in  these  combustion  models  using  a  standard  pre¬ 
sumed  probability  distribution  function  (PDF)  approach  [10,  14,  19].  This  approach  assumes 
that  the  sub-filter  distribution  of  the  independent  flamelet  parameters  may  be  modeled  us¬ 
ing  presumed  distributions  whose  exact  shapes  depend  on  the  LES  filtered  moments  of  the 
parameters.  The  filtered  value  of  any  chemical  quantity  may  then  by  determined  by  inte¬ 
grating  the  product  of  the  flamelet  solution  and  the  presumed  distribution.  In  both  the 
strained  and  unstrained  models,  a  beta-PDF  (/?)  will  be  used  to  describe  sub-filter  progress 
variable  distributions.  A  beta-PDF  depends  on  a  mean  and  a  variance,  so  that  filtered 
quantities  in  the  unstrained  model  are  calculated  as 

Me,  c"5)  =  J  mc)  •  P(C-  c,  a12)  dc.  (51) 

Beta-PDFs  typically  describe  distributions  that  are  bounded  between  zero  and  one,  but 
the  maximum  progress  variable  in  the  flame  is  less  than  unity.  To  account  for  this,  the 
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beta-PDF  is  rescaled  so  that  it  describes  distributions  between  C  =  0  and  the  maximum 
expected  progress  variable  value,  C  =  0.20. 

In  the  strained  model,  the  mass  fraction  of  hydrogen  is  employed  as  an  additional 
tabulation  coordinate.  This  mass  fraction  is  used  to  determine  which  strained  flamelet 
solution  from  the  chemistry  database  is  accessed  by  the  flow  solver.  Just  as  in  the  unstrained 
model,  the  mean  and  variance  of  the  progress  variable  in  any  LES  cell  can  be  used  to  specify  a 
beta-PDF,  /3(C]  C,  C"2),  that  approximates  the  distribution  of  C  in  the  cell.  The  beta-PDF 
can  be  integrated  against  the  hydrogen  profile  from  any  single  strained  flamelet  solution  to 
determine  its  conditional  hydrogen  mass  fraction,  Yjj\C.  The  particular  strained  flamelet 
that  is  selected  for  access  is  the  one  whose  conditional  hydrogen  mass  fraction  matches  that 
of  the  transported  LES  hydrogen  mass  fraction,  Yu.  This  comparison  is  mathematically 
performed  by  presuming  that  the  distribution  of  the  conditional  hydrogen  value  is  a  delta- 
PDF.  Chemical  quantities  then  are  written 

MC,  C^,  Yh)  =  JI  UC,  Yh)/3(C ;  C,  cF2)5(YH  -  Yh\C)  dC  dYH.  (52) 

The  variance  of  the  progress  variable  is  calculated  using  an  algebraic  model  with  a  dynam¬ 
ically  computed  coefficient  Cs  [23],  C"2  =  CsA2|VCj2. 

2.5.4  LES  Approach 

LES  computations  of  the  slot-jet  DNS  are 
performed  using  two  different  mesh  resolu¬ 
tions.  The  first  mesh  consists  of  1.2  million 
cells,  and  the  corresponding  filter  width  ra¬ 
tios  are  A/77  =  8  and  A /Ip  =  0.5.  The  second 
LES  mesh  consists  of  9.3  million  cells,  and  the 
corresponding  filter  ratios  are  A/77  =  4  and 
A /Ip  =  0.25.  Although  they  are  relatively 
fine,  these  LES  resolutions  rigorously  test  the 
models  because  they  describe  unresolved  in¬ 
ner  premixed  flame  structures  using  only  the 
flamelet  databases. 

Both  meshes  employ  stretching,  so  that 
the  cell  size  increases  in  the  downstream 
and  radial  directions.  The  LES  simulations 
are  performed  using  a  parallel  finite  differ¬ 
ence  code  [56]  that  is  advanced  in  time  using 
a  Crank-Nicolson-type  second  order  implicit 
scheme.  Spatial  gradients  are  calculated  us¬ 
ing  second  order  schemes  for  velocities  and 
third  order  schemes  for  scalars.  The  code  is 


Figure  27:  Schematics  of  the  strained  model  LES 
runs.  The  C  =  0.065  surface  is  plotted  over  contours 
of  the  transported  Yh  scalar.  Left:  1.2  million  cells. 
Right:  9.3  million  cells. 
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run  in  a  low-Mach  formulation,  and  a  Poisson  equation  is  solved  for  the  pressure  variable 
that  enforces  mass  conservation. 

Both  the  strained  and  unstrained  flamelet  models  are  applied  on  the  1.2  million  cell 
mesh,  while  only  the  strained  model  is  applied  on  the  9.3  million  cell  mesh.  A  schematic  of 
the  strained  flamelet  results  on  the  two  different  meshes  is  shown  in  Fig.  27.  Comparison 
of  the  flame  fronts  in  Figs.  24  and  27  indicates  that  the  features  of  the  flame  surface  are 
reasonably  captured  by  the  LES.  The  hydrogen  mass  fraction  that  is  shown  in  the  contour 
plots  in  Fig.  27  on  a  streamwise-transverse  plane  is  characterized  by  relatively  smooth 
structures  on  the  burned  side  of  the  front  due  to  the  particularly  high  diffusivity  of  hydrogen 
at  high  temperatures. 

C  and  Yh  scalars  are  both  solved  for  in  the  LES  using  standard  transport  equations. 
The  Lewis  number  of  the  progress  variable  is  set  as  Lee  =  1;  while  that  of  the  hydrogen 
radical  is  set  as  Le#  =  0.18. 

It  is  known  [49]  that  premixed  flame  fronts  are  under-resolved  in  standard  implicitly 
filtered  LES  computations.  This  issue  is  dealt  with  in  the  current  study  by  coupling  the 
progress  variable  transport  equation  with  a  level  set  solution  according  to  the  procedure  in 
reference  [14].  In  this  coupling  the  level  set  solution  is  used  to  adjust  the  progress  variable 
chemical  source  in  the  immediate  vicinity  of  the  flame  front.  This  adjustment  ensures  that 
the  flame  front  propagates  at  an  appropriate  turbulent  burning  velocity  even  in  the  presence 
of  numerical  errors  caused  by  poor  flame  resolution.  The  source  term  in  the  hydrogen  mass 
fraction  transport  equation  is  not  directly  coupled  with  the  level  set  because  it  rapidly 
responds  to  the  progress  variable  field  and  the  changes  that  the  level  set  induces  in  this 
field. 

The  transported  level  set  requires  a  laminar  burning  velocity,  a  model  for  the  turbulent 
burning  velocity,  and  a  method  of  reinitializing  the  level  set  field  to  ensure  it  maintains 
a  smooth  gradient.  The  laminar  burning  velocity  is  accessed  from  the  unstrained  and 
strained  flamelet  databases  using  Eq.  (51)  and  Eq.  (52),  respectively.  The  influence  of 
the  local  hydrogen  mass  fraction  is  therefore  accounted  for  in  the  strained  model  LES 
simulations.  The  model  of  Pitsch  [49]  is  used  to  describe  the  influence  of  turbulence  on  the 
burning  velocity,  and  level  set  reinitialization  is  accomplished  using  a  parallel  fast  marching 
method  [57]. 

2.5.5  Results  And  Discussion 

Progress  Variable 

Time  averaged  progress  variable  and  progress  variable  source  term  results  from  the 
strained  and  unstrained  LES  runs  are  compared  with  corresponding  DNS  results  in  Fig.  28. 
At  the  first  three  measurement  stations,  the  unstrained  flamelet  LES  (black  line)  is  seen  to 
over-predict  the  progress  variable  source  term.  The  stretched  flamelet  models,  conversely, 
are  able  to  predict  how  turbulent  mixing  perturbs  the  flame  structure  and  lowers  the  source 
terms.  These  different  conditional  source  term  profiles  do  not  affect  the  progress  variable 
solutions  near  the  nozzle,  where  resolved  mixing  processes  control  scalar  transport.  The 
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Figure  28:  Progress  variable,  progress  variable  source  term,  and  unconditional  and  conditional  hydrogen 

from:  DNS  (•  •  •);  unstrained  flamelet  LES  ( - );  strained  flamelet  LES  with  1.2  million  cells  ( - ); 

strained  flamelet  LES  with  9.3  million  cells  ( - ). 


influence  of  the  source  term  does  appear,  however,  by  the  time  fluid  parcels  reach  the 
X/H=7.0  station.  There  the  unstrained  model  over-predicts  the  progress  variable,  indicating 
that  the  burning  velocity  is  not  correctly  described.  This  error  in  the  C  predictions  persists 
at  the  farthest  downstream  station  (X/H=ll).  In  contrast  with  the  unstrained  model,  the 
strained  model  is  able  to  reproduce  both  the  source  term  and  progress  variable  DNS  results. 
Some  mesh  sensitivity  is  noticeable  at  X/H=ll,  but  agreement  with  the  DNS  improves  with 
increasing  LES  resolution,  as  expected. 

Hydrogen 

LES  and  DNS  hydrogen  results  are  compared  in  Fig.  28,  where  both  unconditioned  and 
conditional  data  are  plotted.  The  unstrained  model  tends  to  overpredict  Yu  near  the  inlet, 
where  turbulence  induced  strain  is  strongest.  The  strained  flamelet  models,  conversely, 
agree  well  with  the  DNS.  This  trend  is  consistent  with  the  progress  variable  data  in  Fig.  28 
in  that  strain  tends  to  moderate  the  progress  variable  source  term  and  the  buildup  of  the 
radical  pools  that  drive  flame  propagation.  At  downstream  locations  where  turbulence 
has  attenuated,  Fig.  28  indicates  that  the  unstrained  model  predicts  conditional  hydrogen 
profiles  quite  well.  The  unconditional  data  is  poorly  predicted,  however,  due  to  the  over¬ 
prediction  of  C  that  was  discussed  above. 

Minor  Species 

Minor  species  data  from  the  LES  and  DNS  runs  are  compared  in  Fig.  29.  The  strained 
model  has  more  difficulty  predicting  CO  than  it  does  predicting  the  transported  C  and  Yu 
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Figure  29:  CO  (left)  and  OH  (right)  mass  fractions  from:  DNS  (•  •  •);  unstrained  flamelet  LES  ( - ); 

strained  flamelet  LES  with  1.2  million  cells  ( - );  strained  flamelet  LES  with  9.3  million  cells  ( - ). 


scalars.  The  first  two  stations  in  Fig.  29,  for  example,  are  characterized  by  significant  CO 
over-predictions.  Thornber  et  al.  [58]  also  observed  some  error  in  their  LES-CMC  predictions 
of  CO  near  the  X/H=2.0  station,  although  the  error  changed  from  an  over-prediction  of  CO 
on  a  coarse  mesh  to  an  under-prediction  on  a  finer  mesh.  These  errors  may  arise  because 
CO  is  responsive  to  the  details  of  the  imposed  strain  field.  A  reduced  flamelet  manifold  in 
which  this  species  is  assumed  to  be  completely  described  by  two  parameters  may  therefore 
be  insufficient.  The  progress  variable,  which  contains  CO,  benefits  from  the  fact  that  it  is 
defined  to  simultaneously  consider  CO  and  CO2,  reducing  sensitivity  to  oxidation  reactions 
involving  both  species. 

At  downstream  stations  the  quality  of  the  strained  model  CO  predictions  improve,  and 
little  sensitivity  to  the  LES  mesh  is  noticeable.  Strained  flamelet  predictions  of  the  OH 
radical  are  in  relatively  good  agreement  with  the  DNS  throughout  the  flame,  while  the 
unstrained  model  OH  predictions  are  subject  to  errors  at  downstream  locations. 

2.5.6  Summary 

This  study  proposed  a  strained  flamelet  model  for  describing  high  Karlovitz  number  pre¬ 
mixed  flames  in  the  context  of  LES.  The  model  was  validated  by  performing  a  series  of 
LES  computations  of  a  premixed  slot-jet  flame  DNS  that  exists  on  the  edge  of  the  broken 
reaction  zones  regime.  The  LES  computations  compared  the  strained  flamelet  model  with 
an  unstrained  model,  and  considered  two  LES  mesh  resolutions.  It  was  found  that  strain 
effects  had  a  leading  order  influence  on  flame  propagation,  and  that  the  strained  flamelet 
model  could  accurately  predict  this  influence  in  LES. 
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3  Chemical  Mechanisms  For  Jet  Fuel  Surrogates 

3.1  Review  Of  Objectives 

Part  two  of  this  project  aimed  to  develop  a  consistent,  reliable,  and  compact  chemical 
scheme  for  a  potential  jet  fuel  surrogate.  The  first  project  objective  was  the  formulation  of 
a  surrogate  mixture  that  can  accurately  represent  the  average  properties  of  jet  fuel.  This 
formulation  is  accomplished  using  a  constrained  optimization  approach,  as  discussed  below. 
The  second  objective  was  the  development  of  a  chemical  mechanism  that  can  describe  the 
oxidation  of  this  potential  surrogate.  Development  of  this  mechanism  is  discussed  in  the 
second  half  of  the  section. 

A  single  consistent  chemical  mechanism  was  developed  earlier  [59]  to  describe  the  oxi¬ 
dation  of  smaller  hydrocarbons,  from  Ci  to  Cs  species.  Pathways  for  soot  formation  were 
included  in  its  formulation.  This  mechanism  is  used  as  a  base  model  for  the  kinetic  approach 
that  is  developed  here.  In  the  current  project,  the  base  model  is  extended  so  that  it  includes 
the  oxidation  pathways  of  the  individual  components  of  the  surrogate  mixture.  The  com¬ 
bined  chemical  scheme’s  description  of  oxidation  is  validated  at  each  step  of  this  process. 
The  resulting  chemical  model’s  ability  to  describe  the  global  combustion  characteristics  of 
jet  fuel  is  then  examined. 

The  mechanism  development  process  utilizes  the  multi-stage  reduction  strategy  devel¬ 
oped  under  a  previous  AFOSR-funded  project.  The  reduction  tools  exploited  include  the 
Directed  Relation  Graph  with  Error  Propagation  (DRGEP)  reduction  method  [60].  DRGEP 
automatically  eliminates  from  a  kinetic  mechanism  the  species  and  reactions  that  do  not  con¬ 
tribute  significantly  to  the  overall  dynamics  of  a  chemical  process.  Additionally,  a  chemical 
lumping  scheme  [61]  is  used  to  group  chemical  isomers  into  lumped  species.  All  numerical 
calculations  have  been  performed  using  the  FlameMaster  program  [62], 

3.2  Definition  Of  Surrogates  For  Jet  Fuels 
3.2.1  Choice  Of  Individual  Components 

A  natural  procedure  to  select  suitable  components  of  the  surrogate  mixture  for  a  particular 
fuel,  is  to  identify  a  representative  hydrocarbon,  to  stand  for  each  of  the  major  hydrocarbon 
classes  found  in  the  real  fuel.  Based  on  molecules  identified  as  relevant  to  jet  fuels  [8],  the 
components  of  the  jet  fuel  surrogate  for  this  work  have  been  chosen  as:  (a)  n-dodecane,  to 
represent  the  paraffin  class;  (b)  methylcyclohexane,  to  represent  the  naphthene  class;  and 
(c)  m-xylene  to  represent  the  aromatics. 

This  choice  is  motivated  by  several  observations.  First,  n-dodecane  has  physical  proper¬ 
ties  close  to  JP-8/Jet-A  over  temperature  ranges  of  100-650°C  [63,  64],  Second,  m-xylene 
aids  the  surrogate  in  reproducing  the  sooting  tendencies  of  a  jet  fuel.  Finally,  methylcy¬ 
clohexane  is  the  simplest  branched  naphthene  that  can  be  modeled  reliably.  The  global 
ignition  and  flame  propagation  characteristics  of  these  components  have  been  examined  in 
several  experimental  studies,  and  some  of  their  key  chemical  reaction  pathways  have  also 
been  the  object  of  theoretical  and  experimental  kinetic  rate  determinations.  It  is  important 
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Parameters 

Gasoline 

Jet  fuel 

Diesel  fuel 

Lower  HV  [MJ/kg] 
Carbon  number  range 
Approximate  formula 
Liquid  density  [kg/1] 
Molecular  weight  [g/mol] 

43.4 

4-12 

C6.9H13.5 

0.735 

~  96.3 

43.3 

8-16 

C11H21 

0.775-0.840 

~  153 

42.7 

9-23 

Ci6H28 

0.850 

~  220 

Table  1:  Average  properties  of  transportation  fuels  key  to  define  surrogates  for  combustion  applications. 
Data  compiled  from  several  sources  [63,  65,  66,  67]. 

to  choose  representative  components  that  have  been  carefully  studied,  in  order  to  assess  the 
ability  of  the  surrogate  kinetic  model  for  the  individual  component  description,  which  is  in 
itself  key  to  the  performance  of  the  multi-component  surrogate  model. 

3.2.2  Composition  Of  Individual  Components 

A  constrained  optimization  approach  is  used  to  determine  the  composition  of  the  individ¬ 
ual  components  of  the  surrogate  mixture.  Given  a  choice  of  components  to  make  up  the 
surrogate  (here  n-dodecane,  methylcyclohexane,  and  m-xylene),  an  optimized  component 
composition  is  determined,  so  that  the  properties  of  the  surrogate  fuel  resemble  the  target 
real  fuel  properties.  The  procedure  itself  was  developed  as  a  part  of  a  previous  AFOSR 
grant.  A  brief  description  is  provided  here. 

The  average  properties  of  typical  transportation  fuels  that  need  to  be  considered  to 
design  a  surrogate  for  combustion  applications  are  shown  in  Table  1.  The  problem  of  de¬ 
termining  the  composition  of  individual  components  in  the  surrogate  mixture  is  formulated 
as  a  constrained  optimization  problem.  The  average  real  fuel  properties  are  the  desired  op¬ 
timization  targets,  and  some  of  the  target  real  fuel  properties  could  be  used  as  constraints. 
In  recasting  the  problem,  the  mixture  properties  are  determined  by  exploiting  the  fact  that 
most  of  these  target  real  fuel  properties  are  indeed  bulk  properties;  the  multi-component 
surrogate  fuel’s  average  properties  are  hence  expressed  as  combinations  of  individual  com¬ 
ponent  properties.  Structural  group  analysis  is  also  used  wherever  appropriate,  for  instance 
to  determine  Threshold  Sooting  Index  (TSI)  of  the  surrogate  fuel. 

The  optimal  component  composition  of  a  jet  fuel  surrogate  that  is  obtained  by  solving 
the  constrained  optimization  problem  for  n-dodecane,  methylcyclohexane,  and  m-xylene  is 
provided  in  Table  2.  It  can  be  seen  that  the  proposed  jet  fuel  surrogate  agrees  with  the 
target  real  fuel  properties  as  far  as  the  composition  of  major  hydrorcabons,  and  H/C  ratio 
(which  is  indicative  of  the  heating  value  of  the  jet  fuel)  is  concerned.  The  sooting  tendency 
measure,  Threshold  Sooting  Index ,  of  the  jet  fuel  is  also  captured  by  the  model  surrogate. 
However,  the  average  chemical  formula  and  the  fuel  molecular  weight  are  different  between 
the  real  fuel  and  the  model  surrogate.  If  a  heavier  naphthene  were  to  be  included,  these 
differences  could  be  resolved.  In  spite  of  these  discrepancies,  it  will  be  shown  in  the  later 
sections  that  the  surrogate  proposed  here  accurately  describess  the  real  fuel’s  combustion 
characteristics. 
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Jet-A/JP-8  properties 

Model  JP-8  Surrogate 

H/C  ratio 

1.91  ±  0.05 

1.94 

Average  formula 

C11H21 

C9.64Hi8.75 

Liquid  density  (kg/1) 

0.810 

0.769 

Molecular  weight  (g/mol) 

153 

134.80 

Threshold  Sooting  Index  (TSI) 

15 

15.34 

~60%  paraffins 

62%  n-dodecane 

Composition  (%  volume) 

~20%  cyclo  paraffins 

19.8%  methylcyclohexane 

~18%  aromatics 

17.6%  m-xylene 

Table  2:  A  model  jet  fuel  surrogate 


3.3  A  Single  Chemical  Mechanism  Describing  A  Jet  Fuel  Surrogate 

Having  determined  the  composition  of  an  optimal  jet  fuel  surrogate  as  described  above,  it 
is  our  objective  to  elaborate  the  development  of  a  single,  consistent  and  reliable  chemical 
scheme  to  accurately  model  such  a  surrogate  mixture.  The  basic  idea  is  pictorially  rep¬ 
resented  in  Fig.  30.  As  mentioned  earlier,  the  present  work  uses  the  detailed  mechanism 
published  by  Blanquart  et  al.  [68]  as  the  base  chemical  model.  This  mechanism  has  been 
developed  for  the  oxidation  of  thirteen  fuels  ranging  from  the  Ci  to  the  Cs  species,  and 
including  alkanes  such  as  n-heptane  and  iso-octane.  Additionally,  aromatic  species  such 
as  benzene  and  toluene  are  included.  This  base  mechanism  has  been  extensively  validated 
against  ignition  delay  times  and  laminar  burning  velocities  over  a  wide  range  of  temper¬ 
atures  and  pressures.  Finally,  it  has  been  applied  to  a  series  of  laminar  premixed  and 
diffusion  flames.  In  all  the  cases  investigated,  the  mechanism  was  found  to  reproduce  the 
experimental  measurements  well. 

Here,  a  detailed  discussion  of  the  extension  of  the  base  chemical  model  is  provided.  This 
extension  will  introduce  the  reaction  pathways  of  several  substituted  aromatics,  including 
m-xylene.  Particular  attention  is  paid  to  demonstrate  the  ability  of  the  mechanism  to  pre¬ 
dict  the  kinetics  of  m-xylene,  since  it  is  a  component  of  the  proposed  surrogate.  Thereafter, 
the  procedure  to  extend  the  chemical  model  to  include  kinetics  of  n-dodecane  is  discussed, 
and  extensive  validation  tests  for  this  normal  alkane  are  provided.  A  similar  approach  to 
extend  the  chemical  scheme  to  include  reaction  pathways  of  the  naphthene  representative  in 
the  proposed  surrogate,  methylcyclohexane,  is  then  discussed.  The  resulting  chemical  mech¬ 
anism  would  serve  to  reliably  predict  the  kinetics  of  the  model  jet  fuel  surrogate  proposed 
in  Table  2. 

3.3.1  Extension  To  Substituted  Aromatics 

Aromatic  compounds  are  major  constituents  of  real  engine  fuels:  25%  by  volume  in  gaso¬ 
line,  33%  in  diesel,  and  16%  in  jet  fuels  (JP-8,  Jet  A/A-l)  [69].  They  are  used  as  anti-knock 
additives  to  enhance  the  octane  number  of  the  fuels,  since  they  have  high  resistance  to 
auto- ignition  [70].  Aromatic  species  also  play  a  crucial  role  in  the  formation  of  soot  as 
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Figure  30:  Schematic  of  the  development  of  a  chemical  model  for  the  oxidation  of  a  proposed  jet  fuel 
surrogate. 


they  enhance  the  formation  of  soot  precursors  such  as  Polycyclic  Aromatic  Hydrocarbons 
(PAHs).  Therefore,  understanding  and  accurately  modeling  the  chemistry  of  aromatic  com¬ 
pounds  is  an  essential  component  of  the  process  of  describing  surrogate  fuels. 

Much  work  in  the  past  has  been  devoted  to  modeling  combustion  of  toluene.  Chemical 
mechanisms  for  high-temperature  (875  -  1500  K)  gas-phase  oxidation  of  benzene,  toluene, 
ethylbenzene,  and  propylbenzene  are  described  and  discussed  in  detail  by  Brezinsky  et 
al.  [70,  71].  A  detailed  chemical  kinetic  mechanism  for  combustion  of  toluene  at  intermediate 
and  high  temperatures  has  been  assembled  and  evaluated  by  Lindstedt  and  Maurice  [72] 
for  a  wide  range  of  oxidation  regimes. 

Mechanisms  for  dimethylbenzenes  (xylenes)  have  also  been  established  in  the  literature. 
Battin-Leclerc  et  al.  [73]  measured  ignition  delay  times  of  xylenes  in  a  shock  tube  and 
proposed  a  reaction  scheme  to  reproduce  the  experimental  results.  Gail  and  Dagaut  [74] 
studied  oxidation  of  p-xylene  (1,4-dimethylbenzene)  and  m-xylene  (1,3-dimethylbenzene) 
in  a  Jet  Stirred  Reactor  (JSR)  at  atmospheric  pressure.  Similarly,  a  detailed  kinetic  model 
was  developed  to  describe  the  oxidation  of  p-xylene  that  represents  the  experimental  JSR 
results.  Finally,  a  detailed  chemical  mechanism  was  proposed  for  the  oxidation  of  o- xylene 
(1,2-dimethylbenzene)  and  m-xylene  by  the  same  group  [75,  76]  to  represent  their  JSR  and 
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ignition  delay  time  results. 

In  the  study  of  bi-ringed  substituted  species,  Shaddix  et  al.  [77]  proposed  reaction  path¬ 
ways  for  the  high  temperature  oxidation  of  1-nrethylnaphthalene  based  on  the  observed 
intermediate  species  profiles  in  plug  flow  reactor  (PFR)  measurements  [77].  Pitsch  [78]  pro¬ 
posed  a  mechanism  for  the  oxidation  of  1-methylnaphthalene  and  validated  the  mechanism 
with  PFR  and  ignition  delay  time  data.  More  recently,  Mati  et  al.  [79]  also  developed  a 
model  for  1-methylnaphthalene  oxidation  and  validated  their  mechanism  with  JSR  mea¬ 
surements  and  ignition  delay  times. 

As  is  evident,  in  almost  all  of  the  work  done  thus  far,  the  focus  has  been  on  developing 
detailed  chemical  mechanisms  for  the  oxidation  of  individual  fuel  species.  In  the  present  ef¬ 
fort,  a  single  mechanism  has  been  developed  from  on  the  base  model  described  earlier.  This 
single  mechanism  can  be  used  to  describe  several  substituted  aromatics:  toluene,  ethylben¬ 
zene,  styrene,  a-nrethylnaphthalene  and  m-xylene  [80].  The  focus  of  this  report  is  laid  on 
m-xylene  since  this  is  a  component  of  the  proposed  surrogate  for  jet  fuels.  The  m-xylene 
validation  tests  that  will  be  discussed  include  ignition  delay  times,  species  concentration 
profiles  in  shock  tube  experiments,  PFR  data,  and  laminar  burning  velocities. 

For  clarity,  in  the  rest  of  the  article,  abbreviations  according  to  those  introduced  by 
Frenklach  et  al.  [81]  have  been  used.  For  instance,  Ai  is  benzene,  A1CH3  is  a  methyl- 
substituted  species  consisting  of  one  aromatic  ring  (i.e.  toluene),  Ai (0113)2  refers  to  a 
bi-methyl  substituted  one-ringed  aromatic  (i.e.  xylene),  Ai-  refers  to  a  phenyl  radical, 
A1CH3  refers  to  a  methylphenyl  radical,  A2-  is  the  naphthyl  radical,  and  A2CH3  refers  to 
a  methyl-substituted  species  consisting  of  two  aromatic  rings  (i.e.  methylnaphthalene) .  In 
the  following  sections,  substitutions  on  a  naphthyl  ring  refer  to  the  a  site  if  not  stated 
otherwise. 

Mechanism  Development 

The  present  chemical  mechanism  is  composed  of  158  species  and  1804  reactions  (forward 
and  backward  reactions  counted  separately) .  Most  of  the  reactions  are  treated  as  reversible 
with  the  rate  constants  for  the  reverse  reactions  evaluated  from  the  corresponding  equilib¬ 
rium  constants. 

The  oxidation  pathways  of  substituted  aromatic  molecules  are  strongly  coupled.  For 
instance,  the  chemistry  of  xylene  depends  on  the  chemistry  of  toluene,  which  in  turn  depends 
on  the  chemistry  of  benzene.  The  present  chemical  model  proposes  a  consistent  framework 
to  describe  the  pyrolysis  and  oxidation  of  these  aromatic  fuels.  Consistency  here  has  three 
different  meanings.  First,  the  chemistry  of  all  the  fuels  is  described  by  one  single  kinetic 
mechanism.  This  means,  for  example,  that  the  toluene  part  of  the  mechanism  is  used  to 
describe  the  oxidation  of  toluene  itself,  but  also  as  a  sub-mechanism  in  the  oxidation  of 
larger  species,  such  as  xylene,  ethylbenzene,  and  others.  In  addition,  the  chemistry  for  all 
aromatic  fuel  molecules  is  based  on  the  same  chemical  mechanism  for  smaller  hydrocarbons. 
Second,  the  rate  constants  used  for  reactions  of  a  given  class  of  chemical  reactions  are 
correlated.  When  published  data  is  unavailable,  rates  have  been  determined  by  analogy  to 
other  reactions.  The  underlying  rates  used  in  these  analogies  are  consistent  throughout  the 
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mechanism.  As  an  example,  the  rate  of  H-abstraction  from  a  methyl  group  in  toluene  and  in 
xylene  are  taken  in  analogy  to  the  same  reaction  and  are  not  independently  adjusted.  Third, 
none  of  the  reaction  rates  has  been  fitted  to  match  experimental  data.  All  kinetic  rate  data 
are  consistent  with  published  literature.  Note  that  while  this  mechanism  is  developed  to 
describe  combustion  of  several  components,  which  might  be  crucial  for  surrogates  of  realistic 
transportation  fuels,  blends  of  fuels  have  not  been  validated  here.  The  term  consistent 
therefore  does  not  refer  to  the  interactions  of  oxidation  reactions  of  different  fuels. 

The  reaction  rate  parameters  used  in  the  mechanism  are  obtained  from  quantum  chem¬ 
ical  calculations  reported  in  the  literature  or  from  experimental  measurements  from  the 
literature  when  available.  In  certain  cases,  the  reaction  rates  are  adapted  from  the  rates  of 
reactions  of  the  smaller  hydrocarbon  species,  which  for  consistency  have  been  taken  from 
the  base  mechanism  described  by  Blanquart  et  al.  [68]. 

Xylene  Oxidation 

The  chemistry  of  dimethylbenzene  or  xylene  (Ai (0113)2)  oxidation  in  this  reaction  scheme 
is  mostly  based  on  the  chemistry  of  toluene.  The  reaction  scheme  has  been  developed  for 
to- xylene  oxidation,  and  the  reaction  pathways  included  are  primarily  those  of  importance 
to  m-xylene.  The  main  reaction  pathways  involving  the  radicals  and  stable  intermediates 
formed  during  the  oxidation  of  xylene  are  shown  in  Fig.  31. 


Figure  31:  Reaction  pathway  for  m-xylene  oxidation. 


Xylene  Decay  And  Xylyl  Decomposition 

Following  Erndee  et  al.  [82] ,  m-xylene  is  oxidized  by  sequential  reaction  and  removal  of  the 
methyl  side  chains.  Abstraction  of  a  benzylic  H  from  xylene  forms  the  xylyl  (methylbenzyl) 


59 


o* 

Figure  32:  O  attack  on  the  aromatic  ring  of  xylene. 


radical  (A1CH3CH2)  and  the  substitution  of  a  methyl  group  by  an  H  atom  produces  toluene 
and  a  methyl  radical 


Ai(CH3)2  +  H^  AiCH3CH2  +  H2  (53) 

^  AiCH3  +  CH3  .  (54) 

Substitution  of  a  methyl  group  in  xylene  by  an  OH  radical  leads  to  the  formation  of  cresols 
(HOAiCH3).  Xylene  also  leads  to  the  formation  of  methylphenyl  radicals  upon  loss  of  a 
methyl  group. 

An  O  atom  attack  on  xylene  leads  to  the  formation  of  a  dimethylphenoxy  radical 
(OAi(CH3)2),  the  bi-substituent  analogue  of  cresoxy.  In  the  present  work,  OAi(CH3)2 
is  assumed  to  form  toluene  with  the  loss  of  CO  and  H  atoms  (see  Fig.  32).  However,  we 
have  not  included  the  OAi(CH3)2  species  and  the  intermediates  in  the  pathway.  Instead, 
we  have  laid  out  a  direct  pathway  to  form  toluene  from  an  O  attack  on  xylene  assigning  the 
full  entrance  channel  rate  to  the  products.  Bypassing  this  intermediate  oxygenated  species 
is  supported  by  the  work  of  Shaddix  et  al.  [77],  who  observed  for  the  1-methylnaphthalene 
+  O  reaction  that  the  intermediate  OA2CH3  is  short-lived.  The  direct  pathway  laid  out 
to  form  the  products  in  the  above  reactions  is  further  justified  by  results  that  we  obtained 
by  simulations  with  a  more  complete  representation  of  this  pathway.  No  differences  were 
observed  in  the  results  when  the  intermediate  oxygenated  species  shown  in  the  xylene  +  O 
pathway  (OAi(CH3)2)  was  introduced  in  the  mechanism,  as 

Ai(CH3)2  +  O^OAi(CH3)2  +  H  (55) 

OA1(CH3)2^  AiCH3  +  CO  +  H,  (56) 

where  the  reaction  rates  for  these  reactions  were  adapted  from  those  of  the  reactions  of 
benzene  +  O  and  AiO  forming  C5H5  +  CO. 

The  xylyl  radical  undergoes  reactions  similar  to  those  of  the  benzyl  radical.  The  pricipal 
decomposition  pathway  is  shown  in  Fig.  33.  Unimolecular  decomposition  of  the  xylyl  radical 
results  in  the  formation  of  methylcyclopentadienyl  radical  (CsH4CH3)  and  acetylene  (C2H2) 

AiCH3CH2  ->  C5H4CH3  +  C2H2.  (57) 

The  methylcyclopentadienyl  radical  (C.5H4CH3)  formed  here  is  analogous  to  the  cycylopen- 
tadienyl  radical  (C5H5)  formed  in  benzyl  radical  decomposition.  Emdee  et  al.  [82]  proposed 
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Figure  33:  Decomposition  of  xylyl  radical. 


Reactions 

A 

n 

E 

A1CH3CH2  +  0  - 

>  AiCHgCHO  +  H 

4.37-  1018 

-1.34 

6.66 

AiCHgCHs  +  0  - 

>  A1CH3  +  ch2o 

5.99  •  1023 

-2.47 

67.75 

AiCHsCHs  +  0  - 

>  A1CH3  +  HCO 

1.97-  1022 

-2.36 

34.11 

A1CH3CH2  +  02 

-A  AiCHaCHO  +  OH 

1.38-  102 

+2.42 

31.13 

A1CH3CH2  +  02  -A  OA1CH3  +  ch2o 

6.57-  103 

+1.87 

20.93 

Table  3:  Reactions  whose  rates  have  been  derived.  Refer  to  text  for  details.  Rate  coefficients  in  Arrhenius 
form  (k  =  ATnexp(— E/RT));  Units  are  cm3,  K,  mol,  s  and  kj. 

a  fast  rearrangement  of  the  methylcyclopentadienyl  radical  to  give  benzene  and  H  atom. 

C5H4CH3  -A  C6H6  +  H  (58) 

A  global  reaction  has  been  proposed  in  the  current  mechanism  for  the  unimolecular  de¬ 
composition  of  the  xylyl  radical  forming  benzene,  acetylene,  and  H  atom.  Such  a  global 
reaction  has  also  been  used  in  the  work  of  Battin-Leclerc  et  al.  [73].  However,  for  consis¬ 
tency  with  the  toluene  chemistry,  the  rate  assigned  to  xylyl  radical  decomposition  is  that 
of  the  unimolecular  decomposition  of  benzyl  radicals  [83]  used  in  the  present  mechanism. 

Successive  Oxidation  Of  The  Methyl  Groups 

The  xylyl  radical  reacts  with  an  O  atom  to  form  the  A1CH3CH2O  species  which  decom¬ 
poses  to  give  tolualdehyde,  toluene,  and  methylphenyl  radicals  in  a  manner  similar  to  the 
reactions  of  benzoxyl  radicals  (A1CH2O)  in  toluene  oxidation 


A1CH3CH2  +  0  ±5  A1CH3CHO  +  H 

(59) 

A1CH3  +  HCO 

(60) 

±5  A1CH3  +  ch2o  . 

(61) 

The  species  A1CH3CH2O  is  assumed  to  rapidly  decompose  into  products  and  has  been 
bypassed  in  the  A1CH3CH2  +  O  reaction.  The  reaction  rates  for  the  overall  reactions  have 
been  derived  using  the  corresponding  reaction  rates  of  the  O  attack  on  benzyl  radical  and  the 
decomposition  of  the  benzoxyl  radical  (A1CH2O),  assuming  that  the  A1CH3CH2O  species 
exists  in  quasi  steady  state.  The  reaction  rates  derived  here  are  presented  in  Table  3. 
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The  reaction  of  molecular  oxygen  with  a  xylyl  radical  forms  a  methylbenzylperoxy  rad¬ 
ical  which  further  decomposes  into  tolualdehyde  (A1CH3CHO)  or  a  cresoxy  species  in  a 
manner  similar  to  the  reactions  of  the  benzyl  radical 

A1CH3CH2  +  02  ^  A1CH3CHO  +  OH  (62) 

U  OA1CH3  +  CH2O.  (63) 

In  determining  the  rates  for  this  reaction,  an  approach  similar  to  that  for  the  benzyl  + 
O2  reaction  is  used.  Canneaux  et  al.  [84]  computed  the  rate  constant  for  the  four-centered 
isomerization  of  the  methylbenzylperoxy  radical  to  form  the  methyl(hydroperoxy)benzyl 
radical  using  an  elaborated  CASPT2  method.  This  isomerization  reaction  is  the  rate  deter¬ 
mining  step  for  the  formation  of  tolualdehyde  from  methylbenzylperoxy  radicals.  Murakami 
et  al.  [85]  performed  a  study  at  the  CBS-QB3  level  of  theory  to  compute  the  rate  constants 
and  product  branching  ratios  for  o-xylyl  +  O2  and  its  subsequent  reactions. 

The  individual  reaction  rates  to  directly  form  tolualdehyde  and  cresoxy  species  from 
xylyl  +  O2  are  obtained  assuming  the  methylbenzylperoxy  radical  to  exist  in  quasi-steady 
state.  In  this  calculation,  the  forward  and  the  backward  reaction  rates  of  xylyl+  O2  forming 
methylbenzylperoxy  radical,  and  the  reaction  rate  to  form  cresoxy  species  from  methylben¬ 
zylperoxy  radical  are  obtained  from  Murakami  et  al.  [85],  assuming  the  rates  are  roughly 
the  same  for  both  m- xylene  and  0- xylene;  the  reaction  rate  for  the  conversion  of  methylben¬ 
zylperoxy  to  tolualdehyde  is  derived  from  Canneaux  et  al.  [84],  The  final  rate  constant  for 
the  formation  of  tolualdehyde  was  found  to  be  slightly  smaller  than  that  for  benzyl  (about 
30%  lower).  Given  the  sensitivity  of  the  results  to  this  rate  constant,  it  was  found  very 
important  to  evaluate  this  reaction  rate. 

The  reactions  of  tolualdehyde  are  based  on  the  reactions  of  toluene  and  benzaldehyde. 
The  methyl  side  chain  of  tolualdehyde  reacts  like  the  methyl  group  on  toluene  to  give 
methoxybenzyl  (A1CHOCH2)  radicals.  The  formyl  (-CHO)  side  chain  of  tolualdehyde  reacts 
like  the  -CHO  group  on  benzaldehyde  to  give  methylphenyl  radicals.  Tolualdehyde  can  lead 
also  to  the  formation  of  toluene,  benzaldehyde,  cresol,  and  phenoxy  species  via  substitution 


of  the  side  chain  groups  (-CH3,  -CHO)  by  H  and  OH  radicals 

A1CH3CHO  +  H  ±5  A1CH3  +  HCO  (64) 

±5  AiCHO  +  CH3  (65) 

AiCH3CHO  +  OH  ±5  HOA1CH3  +  HCO  (66) 

->  AiO  +  H  +  CO  +  CH3  .  (67) 


The  methoxybenzyl  radical  (A1CHOCH2)  formed  from  tolualdehyde  undergoes  reac¬ 
tions  similar  to  the  benzyl  radical.  The  primary  product  of  the  attack  of  an  O  atom  on 
the  methoxybenzyl  radical  is  phthalaldehyde  (Ai(CHO)2).  Other  pathways  for  the  same 
reaction  lead  to  the  formation  of  benzaldehyde  and  phenyl  radicals  analogous  to  the  O  atom 
attack  on  benzyl  radicals. 

Finally,  pthalaldehyde  is  consumed  by  reactions  similar  to  those  of  benzaldehyde.  Fol¬ 
lowing  the  abstraction  of  an  H  atom,  pthalaldehyde  forms  the  AiCHOCO  species,  which  fur- 
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Ignition  delay  times 

Plug  flow  reactors 

Burning  velocity 

Battin-Leclerc  et  al.  [73] 
Gail  et  al.  [86,  75] 

Shen  et  al.  [87] 

Emdee  et  al.  [82] 

Johnston  and  Farrel  [88] 
Egofopoulous  et  al.  [89] 

Table  4:  Validation  cases  for  m-xylene  undertaken  in  the  present  study. 


ther  forms  benzaldehyde  with  the  loss  of  CO.  Substitution  of  the  -CHO  group  on  pthalalde- 
hyde  by  H  forms  benzaldehyde  and  substitution  by  OH  is  assumed  to  form  phenoxy  radicals. 
With  these  reaction  pathways,  the  oxidation  of  xylene,  hereafter  reduces  to  the  oxidation 
of  toluene  and  its  derivative  species  which  has  been  discussed  in  detail  in  Narayanaswamy 
et  al.  [80]. 

Validation  Tests  For  m-Xylene 

Table  4  lists  the  validation  cases  considered  for  m-xylene  in  the  present  study.  The 
different  experimental  setups  used  for  the  validation  of  the  present  mechanism  were  chosen 
based  on  their  relevance  to  the  modes  of  combustion  typically  found  in  realistic  engines. 

Experimental  data  for  the  oxidation  of  aromatic  fuels  under  engine  conditions  (high 
pressure)  are  limited.  As  a  result,  most  of  the  validation  cases  have  been  performed  at 
atmospheric  pressure.  While  the  chemical  model  has  been  validated  extensively  and  under 
a  large  range  of  conditions,  there  remain  some  uncertainties  on  the  behavior  of  the  model 
in  more  realistic  conditions. 

Ignition  Delay  Results 

Comparison  of  modeled  results  with  the  experimentally  measured  ignition  delay  times 
for  xylene  isomers  are  presented  next.  In  all  our  computations,  the  ignition  delay  times  are 
defined  in  the  same  way  as  in  the  respective  experimental  studies. 

Battin-Leclerc  et  al.  [73]  performed  shock  tube  ignition  measurements  for  all  three  iso¬ 
mers  of  xylene  at  reflected  shock  conditions  of  1330  to  1800  K,  6.7  to  9  bar,  4>  =  0.5  to 
2.0,  and  92  to  98%  argon  dilution.  They  note  that  the  three  isomers  show  almost  iden¬ 
tical  reactivity  at  the  experimental  conditions  investigated.  The  ignition  delay  time  was 
defined  based  on  the  10%  rise  in  the  OH  radical  emission.  Accordingly,  our  computations 
use  the  time  to  reach  10%  of  the  maximum  OH  concentration  to  indicate  the  ignition  delay 
time.  The  comparison  between  the  computed  ignition  times  and  the  experimental  data  for 
m-xylene  is  shown  in  Fig.  34(a). 

As  shown,  the  computed  ignition  delay  times  follow  the  experimental  measurements  at 
cj)  =  2.0  and  those  at  (j>  =  1.0  except  at  very  high  temperatures,  while  they  are  slower  than 
the  experiments  at  lean  conditions.  To  better  understand  the  results  and  properly  estimate 
the  performance  of  the  present  chemical  model,  we  found  it  beneficial  to  consider  also  the 
auto-ignition  behavior  of  the  other  xylenes  isomers. 

Gail  et  al.  [75]  measured  ignition  delay  times  of  o-xylene  at  similar  temperatures  and  fuel 
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(b) 
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Figure  34:  Ignition  delay  times  of  xylenes.  Symbols  -  experiments;  (a)  m-xylene:  Battin-Leclerc  et  al.  [73]; 
(b)  o- xylene:  Gail  et  al.  [75];  (c)  p-xylene:  Gail  et  al.  [86];  (d)  m-xylene:  Shen  et  al.  [87];  solid  lines  -  results 
from  numerical  simulations. 


concentrations  as  investigated  by  Battin-Leclerc  et  al.  [73]  but  at  close  to  atmospheric  pres¬ 
sure.  The  ignition  delay  time  definition  in  these  experiments  and  the  present  computations 
are  based  on  the  time  for  the  maximum  rate  of  rise  of  CH  radical  emission.  Figure  34(b) 
shows  the  comparison  between  the  present  computed  ignition  delay  times  for  m-xylene  and 
the  experimental  measurements.  Note  that  these  are  two  different  isomers  and  that  the 
present  comparison  is  only  motivated  by  the  observation  that  all  the  three  isomers  of  xy¬ 
lene  show  similar  reactivity  at  these  high  temperatures  [73].  From  the  figure,  it  is  seen  that 
the  ignition  delay  times  predicted  by  the  mechanism  are  faster  than  the  experiments,  in 
contrast  to  the  comparison  with  the  data  of  Battin-Leclerc  et  al.  [73]. 

Much  closer  agreement  at  equivalence  ratios  of  1.0  and  1.5  is  seen  in  Fig.  34(c)  where  the 
computed  ignition  delay  times  for  m-xylene  are  compared  with  the  ignition  delay  times  for  p- 
xylene  measured  by  Gail  et  al.  [86]  at  similar  temperature  ranges  and  pressures  spanning  1  to 
1.4  atm.  In  effect,  the  computed  ignition  delay  time  plots  show  approximately  the  same  slope 
for  the  different  experimental  data  sets  at  the  studied  high  temperatures  and  stoichiometric 
conditions.  However,  the  magnitudes  of  the  ignition  delay  shows  some  differences. 
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Finally,  Shen  et  al.  [87]  measured  ignition  delay  times  for  the  three  isomers  of  xylene 
at  lower  temperatures  (941-1408  K)  and  higher  pressures  (9-45  atm)  at  equivalence  ratios 
of  0.5  and  1.0.  Figure  34(d)  shows  the  comparison  between  the  computed  ignition  delay 
times  and  the  experimental  data.  The  ignition  delay  times  are  slightly  overpredicted  by  the 
present  computations  at  both  10  bar  and  35  bar  pressure,  though  the  slope  of  the  variation 
is  predicted  closely.  However,  the  predictions  are  much  closer  to  the  experiments  at  10  bar 
when  compared  to  those  at  the  higher  pressure  of  35  bar.  This  trend  is  consistent  with 
that  observed  for  toluene  ignition  at  elevated  pressures.  This  could  be  expected  because 
the  xylene  oxidation  is  mostly  based  on  the  toluene  chemistry. 

From  the  extensive  comparison  of  ignition  delay  times  for  all  three  isomers  at  various 
temperatures,  pressures,  and  equivalence  ratios,  it  is  possible  to  draw  several  conclusions 
about  the  performance  of  the  chemical  model.  The  mechanism  seems  to  predict  consistently 
longer  ignition  delay  times  than  measured  for  m- xylene  at  elevated  pressures  (from  8  bar 
to  35  bar).  On  the  other  hand,  at  atmospheric  pressure,  the  predicted  ignition  delay  times 
are  shorter  than  measured  experimentally  for  o-xylene  and  p-xylene.  These  results  are 
surprising,  for  the  ignition  delay  times  of  o-xylene  are  expected  to  be  similar  or  shorter 
than  those  of  m-xylene  under  these  conditions.  In  the  absence  of  data  for  m-xylene  ignition 
at  atmospheric  pressure,  it  is  difficult  to  explain  unequivocally  the  observed  discrepancies. 

Results  For  Plug  Flow  Reactor  Experiments 

Concentration  profiles  of  the  fuels  and  their  derivatives  are  validated  against  PFR  exper¬ 
imental  data  next.  The  time  shifts  applied  to  the  experimental  data  (in  order  to  define  zero 
time)  have  been  indicated  in  the  plots. 

Erndee  et  al.  [82]  conducted  atmospheric  flow  reactor  tests  with  m-xylene  and  p- xylene 
at  temperatures  ranging  from  1093  K  to  1199  K  and  equivalence  ratios  from  0.47  to  1.7. 
They  estimate  an  uncertainty  in  mole  fractions  of  ±  5%  for  non-aromatic  carbon-containing 
species  and  ±  10%  for  aromatics  in  their  data.  The  uncertainty  in  H2  is  expected  to  be  ± 
100  ppm  and  ±  5%  for  O2  measurements.  An  absolute  uncertainty  of  ±  10  K  is  expected 
in  measured  temperature. 

The  PFR  simulation  results  obtained  for  m-xylene  from  the  present  scheme  are  validated 
against  the  set  of  experimental  data  for  m-xylene  available  from  Erndee  et  al.  [82] .  Figure  35 
shows  the  comparison  at  stoichiometric  conditions,  1155  K,  and  1  atm. 

The  fuel  (m-xylene)  decay  follows  the  experimental  measurements  closely.  The  predic¬ 
tions  remain  well  within  the  experimental  uncertainty  throughout.  The  decay  pathways  of 
m-xylene  are  similar  to  those  of  toluene.  The  primary  path  of  consumption  is  via  H  atom 
abstraction  from  the  methyl  side  chain  of  m-xylene  forming  methylbenzyl  radicals.  Sub¬ 
stitution  of  a  methyl  group  in  m-xylene  by  an  H  atom  leading  to  the  formation  of  toluene 
contributes  roughly  15%  to  the  fuel  decay.  However,  from  a  reaction  flux  analysis,  this  reac¬ 
tion  is  found  to  be  primarily  responsible  for  all  the  toluene  found  during  m-xylene  oxidation. 
Figure  35a  shows  the  the  computed  profile  of  toluene  in  comparison  with  the  experiments. 
Toluene  leads  to  the  formation  of  benzyl  radicals,  which  on  combining  with  methyl  radicals 
form  ethylbenzene.  The  computed  profile  of  ethylbenzene  in  m-xylene  oxidation  is  shown 
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in  Fig.  35(b). 

Most  of  the  methylbenzyl  radicals  formed  from  to- xylene  result  in  the  formation  of 
tolualdehyde  via  reactions  described  earlier  in  the  article.  Figure  35(c)  shows  that  the 
concentration  profiles  of  tolualdehyde  and  benzaldehyde  are  in  very  good  agreement  with 
the  experiments.  About  70%  of  tolualdehyde  forms  methylphenyl  (A1CH3)  radicals.  It 
was  found  to  be  highly  important  to  include  the  reaction  pathways  that  involve  formation 
of  methylphenyl  radicals  from  tolualdehyde  in  order  to  describe  the  oxidation  of  m- xylene 
fairly  accurately.  The  methylphenyl  radical  further  forms  cresoxy  species  and  fulvene  and 
eventually  leads  to  the  formation  of  benzene.  Figure  35(d)  shows  the  comparison  between 
the  computed  profile  for  benzene  and  the  experimental  measurements.  The  differences 
observed  in  the  benzene  profile  could  be  attributed  to  some  of  the  global  reactions  used  in 
the  mechanism. 

The  concentration  profiles  of  some  of  the  smaller  hydrocarbons  are  also  compared  with 
the  experiments.  The  C2H2,  CH4,  and  H2  profiles  (Figs.  35(d),  35(e))  agree  fairly  well  with 
the  experiments.  The  profile  of  CO  formation  (Fig.  35(f))  and  oxidizer  decay  (Fig.  35(g)) 
are  also  seen  to  agree  closely  with  the  experimental  measurements. 

Laminar  Flame  Speed  Results 

In  this  section,  zero  stretch  burning  velocity  data  computed  for  m-xylene  are  compared 
with  available  experimental  data. 

The  burning  velocities  of  m-xylene  at  different  temperatures  computed  using  the  present 
oxidation  model  are  compared  with  the  experimental  burning  velocity  measurements  for  m- 
xylene  made  by  Johnston  and  Farrel  [88]  in  Fig.  36(c).  It  can  be  seen  that  the  computed 
burning  velocities  are  higher  than  those  reported  in  the  experiments.  However,  such  a  trend 
was  observed  earlier  when  comparing  the  burning  velocities  of  benzene  [68]  and  toluene 
(Fig.  36(a))  and  ethylbenzene  (Fig.  36(b))  with  the  data  of  Johnston  and  Farrel.  There¬ 
fore,  it  was  found  essential  to  validate  against  other  published  data  for  m-xylene  burning 
velocities. 

More  recently,  Ji  et  al.  [89]  measured  laminar  flame  speeds  and  extinction  rates  of  m- 
xylene  as  a  function  of  equivalence  ratio  at  elevated  temperatures  and  atmospheric  pressure. 
Experiments  were  performed  in  the  counterflow  configuration  and  the  flow  velocities  were 
measured  using  Laser  Doppler  Velocimetry.  The  unstretched  laminar  burning  velocity  is 
determined  by  a  linear  extrapolation  of  the  flame  speed  versus  strain  rate  curve.  The 
uncertainty  in  the  experiments  are  about  2  to  3  cm/sec. 

Figure  36(c)  shows  the  comparison  between  the  burning  velocities  computed  using  the 
present  reaction  mechanism  and  the  experimentally  measured  values.  It  can  be  noted 
that  our  computed  values  are  lower  than  those  reported  experimentally  by  about  2  to  4 
cm/s.  Ji  et  al.  [89]  also  report  the  burning  velocities  of  other  fuels.  Their  non-linearly 
extrapolated  burning  velocities  are  found  to  be  lower  than  the  linearly  extrapolated  values 
at  all  equivalence  ratios  by  about  4-5  cm/sec.  This  could  account  for  the  observed  differences 
in  our  results.  It  would  be  interesting  to  compare  our  computed  results  with  the  non-linearly 
extrapolated  burning  velocities  for  m-xylene. 
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Figure  35:  PFR  data  for  m-xylene  at  1155  K,  symbols  -  experiments:  Emdee  et  al.  [82];  solid  lines  -  numerical 
simulations;  shift:  20  ms. 


In  all,  with  the  extensive  validation  undertaken  in  the  present  study,  the  overall  ability 
of  the  model  to  predict  the  oxidation  of  m-xylene  has  been  found  to  be  very  good  given 
the  fact  that  a  single  and  consistent  chemical  model  is  used  to  describe  the  oxidation  of 
a  large  range  of  species.  Similar  satisfactory  agreement  was  also  obtained  for  the  other 
substituted  aromatics  under  consideration,  namely,  toluene,  ethylbenzene,  styrene,  and  a- 
methylnaphthelene.  More  extensive  validation  tests  for  these  substituted  aromatic  species 
are  presented  in  Narayanaswamy  et  al.  [80]. 
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(a)  Toluene 


(b)  Ethylbenzene 


(c)  m-Xylene 


Figure  36:  Laminar  burning  velocity:  symbols  -  experiments;  solid  lines  -  numerical  results. 


3.3.2  Extension  To  n— Dodecane 

Several  detailed  mechanisms  for  n-dodecane  exist  in  literature.  Dahm  et  al.  [90]  and  later 
Herbinet  et  al.  [91],  developed  a  detailed  mechanism  for  n-dodecane,  and  were  able  to  pre¬ 
dict  the  results  of  pyrolysis  experiments  and  thermal  decomposition  in  a  jet-stirred  reactor. 
Ranzi  et  al.  [92]  proposed  and  validated  a  lumped  mechanism  for  n-alkanes  including  n- 
dodecane  for  low  to  high  temperature  oxidation  in  various  configurations.  You  et  al.  [93] 
proposed  a  detailed  kinetic  model  for  normal  alkanes  up  to  n-dodecane  above  850  K.  They 
validated  their  kinetic  model  against  laminar  flame  speeds,  and  ignition  delays  behind  re¬ 
flected  shock  waves,  with  n-dodecane  being  the  emphasis. 

JetSurF  v2.0  [94]  is  an  ongoing  effort  for  jet  fuel  surrogate  mechanism  development  that 
includes  oxidation  of  n-dodecane  at  moderate  to  high  temperatures. 

Of  particular  interest  to  this  work  is  the  detailed  mechanism  proposed  by  Westbrook  et 
al.  [95],  which  includes  high  and  low  temperature  pathways  to  describe  the  pyrolysis  and 
oxidation  of  several  n-alkanes  up  to  n-tetradecane.  As  this  mechanism  was  built  based 
on  well-established  reaction  classes  originally  developed  for  n-heptane,  it  has  been  chosen 
as  the  starting  point  for  our  approach.  This  scheme  for  n-dodecane  is  first  reduced  to  a 
skeletal  level,  and  then  incorporated  in  the  base  mechanism,  now  extended  to  aromatics, 
as  described  in  the  previous  section.  The  methodology  is  discussed  next,  followed  by  a 
demonstration  of  the  performance  of  the  combined  mechanism  for  different  targets. 

Mechanism  Development 

The  reference  mechanism  for  n-dodecane  oxidation  from  Westbrook  et  al.  [95],  referred 
below  as  the  LLNL  mechanism,  has  16313  forward  and  reverse  reactions  among  2112  species. 
First,  this  extensive  mechanism  is  reduced  to  a  skeletal  level  using  a  multi-stage  reduction 
strategy  put  forth  by  Pepiot-Desjardins  and  Pitsch,  involving  species  and  reaction  elimina¬ 
tion  using  the  DRGEP  approach  [60],  and  chemical  species  lumping  [61].  The  database  used 
to  carry  out  the  reduction  includes  plug  flow  reactor-like  configurations,  and  ignition  delays 
for  low  to  high  temperatures  (600-1500  K),  pressures  from  1  to  40  atm,  and  equivalence 
ratios  from  0.5  to  1.5. 

The  skeletal  mechanism  for  n-dodecane  that  is  obtained  from  the  application  of  the 
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(a)  cf>  =  0.5 


(b)  <f>  =  1.0 


1000/T  [K"1] 


(c)  (j>  =  1.5 

Figure  37:  Ignition  delay  times  of  n-dodecane/air,  comparing  the  skeletal  level  mechanism  with  the  detailed 
mechanism  results;  solid  lines  -  detailed  mechanism:  Westbrook  et  al.  [95]  (LLNL). 

DRGEP  approach  consists  of  862  reactions  among  156  species.  Figure  37  compares  igni¬ 
tion  delay  times  computed  using  the  detailed  and  skeletal  schemes  for  lean,  stoichiometric 
and  rich  n-dodecane/air  mixtures  at  20  atm.  The  skeletal  mechanism  provides  very  good 
predictions  in  the  low  (T<850  K)  and  high  (T>1000  K)  temperature  regions.  The  skeletal 
mechanism  predicts  slightly  longer  ignition  delays  in  the  Negative  Temperature  Coefficient 
(NTC)  region.  However,  considering  the  high  reduction  ratio  that  was  achieved  (~95%)  and 
the  small  magnitudes  of  the  errors,  this  skeletal  scheme  is  treated  as  acceptably  accurate 
and  is  used  in  the  subsequent  mechanism  development  steps. 

The  merging  of  the  n-dodecane  with  the  base  mechanism  (including  aromatics)  is  ac¬ 
complished  using  an  interactive  tool  [96]  that  automatically  identifies  common  species  and 
reactions  from  the  different  mechanisms,  and  incompatibilities  between  the  kinetic  data  sets. 
To  ensure  a  smooth  and  consistent  merging,  (a)  Rate  conflicts  detected  during  the  merging 
were  always  resolved  in  favor  of  the  thoroughly  validated  base  mechanism,  therefore  leav¬ 
ing  this  mechanism  virtually  unchanged,  and  (b)  Repeated/redundant  reaction  pathways, 
described  as  elementary  steps  in  the  base  mechanism  while  global  in  the  n-dodecane  sub¬ 
mechanism,  or  vice  versa,  were  removed  from  the  incremental  n-dodecane  sub-mechanism. 
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Shock  tubes 

Burning  velocity 

Ignition  delays 

Species  profiles 

Davidson  et  al.  [97] 
Vasu  et  al.  [64] 

Shen  et  al.  [98] 
Davidson  et  al.  [99] 

Davidson  et  al.  [99] 
Vasu  et  al.  [64] 

Kumar  and  Sung  [100] 

Ji  et  al.  [101] 

Kelley  et  al.  [102] 

Table  5:  Validation  cases  for  n-dodecane  undertaken  in  the  present  study. 


The  combined  mechanism  thus  retains  the  rates  in  the  base  model. 

The  resulting  combined  mechanism  includes  oxidation  pathways  of  n-dodecane  at  a 
variety  of  temperatures.  These  pathways  provide  additions  to  several  hydrocarbon  oxi¬ 
dation  pathways  already  described  in  the  base  mechanism.  The  incremental  n-dodecane 
sub-mechanism  consists  of  267  reactions  among  78  species,  some  of  which  are  lumped.  A 
few  reaction  rate  changes  have  been  introduced  to  improve  the  performance  of  the  mech¬ 
anism  under  different  idealized  configurations  as  described  next.  The  validation  tests  for 
the  substituted  aromatics  presented  in  our  earlier  work  [80]  have  been  repeated  using  the 
combined  mechanism,  and  found  to  produce  negligible  changes  in  results. 

Validation  Tests  For  n— Dodecane 

The  validation  cases  to  evaluate  the  mechanism  focus  on  oxidation,  and  include  ignition 
delays  spanning  wide  ranges  of  temperatures  and  pressures,  species  concentration  profiles 
measured  in  shock  tubes,  and  burning  velocities  obtained  at  different  pressures.  Fuel  pyrol¬ 
ysis  and  transport  dominated  configurations  are  not  tested.  The  list  of  all  validation  cases 
examined  for  n-dodecane  are  summarized  in  Table  5. 

Ignition  Delay  Results 

Ignition  delays  are  computed  using  a  isochor  homogeneous  reactor  configuration,  using 
the  same  ignition  criterion  as  in  the  experiments. 

n-Do decane/ Air  Mixtures 

Vasu  et  al.  [64]  measured  ignition  delays  of  n-dodecane/air  mixtures  over  727-1422  K, 
15-34  atm,  and  0=0.5,  1.0.  Their  experimental  data  scaled  to  20  atm,  and  the  ignition 
delays  computed  using  the  present  mechanism  (labelled  p.w.)  at  20  atm,  600-1500  K  are 
shown  in  Fig.  38.  Some  improvement  over  the  original  LLNL  mechanism  predictions  is  seen, 
due  to  different  H2/O2,  C1-C4  chemistry  in  the  present  base  model. 

Further,  H-abstraction  from  n-dodecane  by  CH3O2  was  introduced  in  the  combined 
mechanism.  This  abstraction  was  analogous  to  a  similar  iso-octane  reaction  pathway  in  the 
base  model,  with  reaction  rates  adapted  from  Carstensen  et  al.  [103].  This  introduction  led 
to  faster  ignition  at  T  >  900  K.  However,  some  differences  are  observed  in  the  NTC  region. 
A  sensitivity  analysis  reveals  that  the  H-abstraction  rate  from  the  fuel  by  HO2  is  key  in 
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(a)  <f>  =  0.5  (b)  <j>  =  1.0 

Figure  38:  Ignition  delay  times  of  n-dodecane/air  mixtures;  symbols  -  experiments:  Vasu  et  al.  [64]. 

determining  the  ignition  delays  around  1000  K: 

n-Ci2H26  +  H02^Ci2H25  +  H202.  (68) 

Alkene  plus  H02  reactions  were  examined  and  found  to  have  negligible  effect  on  the  igni¬ 
tion  delays,  therefore  ruling  out  some  missing  consumption  pathways  involving  H02  in  the 
mechanism  and  prompting  a  more  thorough  investigation  of  reaction  (68). 

While  a  direct  evaluation  of  this  reaction  rate  for  n-dodecane  is  unavailable  in  the 
literature,  it  can  be  estimated  from  kinetic  evaluations  available  for  smaller  carbon  alkanes. 
The  quantum  chemical  calculations  by  Carstensen  et  al.  [103]  and  Aguilera-Iparraguirre  et 
al.  [104]  provide  the  reaction  rates  for  H-abstraction  by  H02  from  primary  and  secondary 
carbons  in  n-butane,  the  latter  study  using  a  more  accurate  theory  than  the  former.  By 
scaling  these  rates  to  account  for  the  number  of  primary  and  secondary  carbons  in  n- 
dodecane,  as  compared  to  n-butane,  two  possible  reaction  rates  for  reaction  (68)  can  be 
obtained.  These  rates  vary  over  a  factor  of  2  to  7  at  different  temperatures,  with  the 
rate  from  Carstensen  et  al.  being  higher.  The  difference  is  attributed  to  the  different 
levels  of  theory  used  in  these  studies  [104],  For  best  agreement  with  experimental  data,  a 
reaction  rate  which  is  1.7  times  that  suggested  by  Aguilera-Iparraguirre  et  al,  well  within 
the  uncertainty  range,  is  chosen  here. 

The  ignition  delays  obtained  with  the  modified  reaction  (68)  are  added  to  Fig.  38,  and 
are  labeled  as  “p.w.  (with  changes).”  As  expected,  the  agreement  between  the  computed 
lines  and  the  experiments  has  improved  significantly  in  the  NTC  region,  with  minimal 
changes  at  high  and  low  temperatures.  Finally,  the  mechanism  captures  the  low  tempera¬ 
ture  data  (T  <  850  K)  at  0=0.5,  but  not  at  0=1.0.  This  will  further  be  investigated  when 
more  measurements  at  such  low  temperatures  become  available. 

Ignition  delays  for  n-dodecane/air  mixtures  were  also  measured  by  Shen  et  al.  [98]  at 
14  and  40  atm,  900-1300  K,  and  0=0.5,  1.0.  In  Fig.  39,  a  comparison  of  the  ignition  delays 
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between  simulations  and  experiments  shows  some  differences  in  all  cases  for  T  <  1050  K, 
with  an  over-prediction  at  14  atm,  and  under-prediction  at  40  atm. 

A  more  comprehensive  comparison  can  be  obtained  by  scaling  Vasu  et  al.  [64]  20  atm 
data  down  to  14  atm  and  up  to  40  atm  using  the  P-1  scaling  law,  used  by  both  groups 
to  scale  data  to  a  common  pressure.  At  14  atm,  T  >  1000  K  and  0=0.5  (Fig.  39(a)), 
the  previously  observed  discrepancies  between  experiments  and  simulations  can  be  traced 
to  potential  outliers  in  the  Shen  et  al.  data  in  this  region.  At  high  temperatures  and 
pressures,  the  computed  results  are  within  the  experimental  uncertainty  and  variability  of 
the  combined  experimental  data. 

Note  that  the  inverse  pressure  law  used  to  scale  the  experimental  data  is  appropriate 
only  in  the  high  temperature  region  (T  >  1000  K).  The  strongly  pressure  dependent  peroxy 
pathways,  dominant  in  the  NTC  region,  would  call  for  a  stronger  relation,  ~P-a,  where 
a>l  [105],  which  would  significantly  improve  the  agreement  between  the  scaled  data  of  Vasu 
et  al.  and  the  computations. 
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Figure  39:  Ignition  delays  of  n-dodecane/air  mixtures;  symbols  -  experiments:  Shen  et  al.  [98]  at  P=14  atm 
and  40  atm  (blue  circles  and  red  squares  respectively);  Vasu  et  al.  [64]  data  at  20  atm  is  scaled  to  14  atm 
and  40  atm  (cyan  and  pink  triangles  respectively)  using  P~x  scaling  (see  text). 


n-Dodecane/02/ 'Argon  Mixtures 

Davidson  et  al.  [99]  and  Vasu  et  al.  [64]  measured  low  fuel  concentration  ignition  times  of 
n-dodecane/02 /argon  mixtures  at  0=1.0,  P=2.25  atm,  and  0=0.5,  P=15  atm  respectively, 
for  T  >  1100  K.  Davidson  et  al.  [97]  also  measured  ignition  delays  of  n- dodecane/02 /Argon 
mixtures  in  an  aerosol  shock  tube,  at  0=0.5  and  P=6.7  atm. 

To  evaluate  the  performance  of  the  mechanism  at  those  different  experimental  conditions 
in  a  consistent  fashion,  the  correlation  developed  by  Horning  et  al.  [106]  for  ignition  delay 
times  of  stoichiometric  n-alkane/02/ Argon  mixtures  is  considered: 

n- alkane  (0  =  1.0)  :  r  =  9.4  x  10“12  P“a55  X0a63  C“a5°  exp(46,  550/RT),  (69) 
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(a)  X°fuel  =  400  ppm 

Figure  40:  Low  fuel  concentration  ignition  delays  of  n-dodecane/02 /argon  mixtures  at  0=1.0;  symbols  - 
experiments:  Davidson  et  al.  [99] .  The  experimental  data,  the  least  square  fit  to  the  experimental  data  (red 
dashed  line)  and  the  computed  ignition  delays  are  scaled  according  to  correlation  from  Horning  et  al.  [106]. 

where  ignition  delays  are  in  microseconds,  pressure  is  in  atmospheres,  Xo2  is  the  oxygen 
mole  fraction,  C  is  the  number  of  carbons  in  the  fuel  molecule,  and  the  activation  energy 
is  in  cal/mole.  This  correlation  was  obtained  from  their  experimental  data  of  propane, 
n-butane,  n-heptane,  n-decane,  for  P=l-6  atm,  T  =  1315-1560  K,  and  Xo2  =  2%-20%. 

The  experimental  conditions  of  Davidson  et  al.  [99],  <0=1.0,  P=2.25  atm,  T=1388- 
1660  K,  Xo2=0.75%  fall  roughly  within  this  validity  range.  In  Fig.  40,  the  experimental 
ignition  delays  measured  by  Davidson  et  al.  [99]  and  the  computed  ignition  delays  have 
been  scaled  according  to  equation  (69),  by  assuming  that  it  holds  true  for  n-dodecane  as 
well.  The  computations  agree  well  with  this  correlation,  though  there  are  differences  when 
compared  to  the  scaled  experimental  data.  While  the  definition  of  ignition  criterion  used 
by  Horning  et  al.  (time  to  reach  peak  CH*  concentration)  and  Davidson  et  al.  (time  to 
reach  50%  peak  OH  concentration)  are  different,  this  should  not  introduce  >10%  difference 
in  ignition  delays.  Consequently,  the  computed  ignition  delays  appear  satisfactory  for  the 
experimental  conditions  studied. 

Ignition  delays  of  lean  n-dodecane/02 /Argon  mixtures  are  discussed  next.  Vasu  et 
al.  [64]  measured  ignition  delays  at  low  fuel  concentrations:  Xn-Ci2H26  =  1000  ppm,  750 
ppm,  514  ppm;  Xo2=10%,  P=15  atm,  <0=0.5,  and  T=1250-1400  K.  To  complement  the 
small  number  of  experimental  data  points,  ignition  delays  of  other  n-alkane/02 /argon  mix¬ 
tures  at  similar  conditions,  expected  to  show  similar  ignition  characteristics  [98],  have  been 
included.  The  experimental  data  for  ignition  delays  of  n- heptane/02 /Argon  mixtures  from 
several  groups  [106,  107,  108]  at  different  experimental  conditions,  included  in  the  study  of 
Horning  et  al.  have  been  scaled  to  15  atm  and  Xo2=10%,  using  p~0-55  and  X^'63  from 
correlation  (69)  discussed  above,  which  has  been  found  to  be  true  for  <0=0.5  as  well  [106]. 
Again,  eq.  (69)  is  assumed  to  be  valid  for  these  conditions. 

In  Fig.  41,  a  comparison  of  the  computed  ignition  delays  with  those  measured  by  Vasu 
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Figure  41:  Low  fuel  concentration  ignition  delays  of  n-dodecane/Ch/argon  at  0=0.5,  P=15  atm;  solid  line  - 
present  work;  symbols  -  experimental  data  [64,  106,  107,  108].  Experimental  data  for  n-heptane/02 /argon 
mixtures  are  scaled  to  the  experimental  conditions  using  correlation  from  Horning  et  al.  [106]. 


et  al.  shows  good  agreement  wherever  data  are  available,  except  at  1000/T  =  0.7  K_1, 
Xn_Ci2H26=1000  ppm.  However,  when  appraising  the  computed  results  together  with  the 
scaled  ignition  data  for  n-heptane,  it  can  be  concluded  that  the  slopes  of  the  computed 
ignition  delays  are  consistent  with  the  scaled  data,  for  all  fuel  concentrations. 

Another  set  of  data  for  ignition  delays  of  lean  n-dodecane/02 /Argon  mixtures  measured 
in  an  aerosol  shock  tube  facility,  at  6.7  atm,  0=0. 5,  T=1050-1330  K  comes  from  Davidson 
et  al.  [97].  In  Fig.  42,  a  comparison  of  the  ignition  delays  of  n-dodecane  between  the 
simulations  and  experiments  shows  wide  disagreement.  However,  in  conjunction  with  the 
ignition  delay  data  for  n-heptane  scaled  to  these  experimental  conditions,  a  few  comments 
can  be  made. 

Ignition  delays  for  n-alkanes  are  expected  to  decrease  as  the  size  of  the  alkane  increases. 
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Figure  42:  Ignition  delay  times  of  n-dodecane/Ch /argon  mixtures  at  P=6.7  atm,  0=0.5;  solid  line  -  present 
work;  symbols  -  experimental  data  [97,  106,  107,  108];  dotted  line  -  fit  to  Davidson  et  al.  data.  Experimental 
data  for  n-heptane/02/argon  mixtures  are  scaled  to  the  experimental  conditions  using  correlation  from 
Horning  et  al.  [106]. 

Accordingly,  the  computed  n-dodecane  ignition  delays  are  shorter  than  the  scaled  ignition 
delay  data  of  n-heptane  throughout  the  high  temperature  region,  which  is  not  the  case  for 
the  experimental  n-dodecane  ignition  data.  Secondly,  the  slope  of  the  computed  ignition 
delays  is  consistent  with  the  scaled  n-heptane  data,  which  is  not  the  case  for  the  n-dodecane 
data  points  (shown  in  Fig.  42).  From  the  extensive  discussions,  consistency  checks  and  the 
range  of  validation  tests  studied,  the  present  mechanism  predicts  the  ignition  delays  of 
n-dodecane  fairly  satisfactorily. 

Species  Time  History  Results 

Davidson  et  al.  [99]  measured  concentration  profiles  of  n-Ci2H26,  C2H4,  OH,  CO2,  and 
H2O  during  the  oxidation  of  n-dodecane/02/ Argon  mixtures  at  P=2.25  atm,  0=1.0.  In 
Fig.  43,  the  species  profiles  computed  at  these  conditions,  are  shown  to  agree  well  with  the 
experiments  for  a  major  part  of  the  studied  time  interval.  Vasu  et  al.  [64]  also  measured 
OH  time  histories  during  the  oxidation  of  n-dodecane/02/Argon  mixtures  at  P=15  atm, 
0=0.5.  In  Fig.  44,  the  computed  OH  profiles  closely  follow  the  measured  rise  of  OH,  except 
at  T=1158  K,  where  the  ignition  seems  to  occur  earlier  than  predicted  by  experiments. 

Laminar  Flame  Speed  Results 

Kumar  and  Sung  [100]  measured  laminar  flame  speeds  of  n-dodecane/air  mixtures  for  a 
range  of  equivalence  ratios,  at  atmospheric  pressure,  and  unburnt  temperature  Tu=400  K 
and  470  K.  Ji  et  al.  [101]  also  measured  flame  speeds  at  P=1  atm  and  Tu=403  K.  Fig¬ 
ure  45(a)  compares  computed  flame  speeds  with  the  experimental  data  of  Kumar  et  al.  and 
the  linearly  and  non-linear ly  extrapolated  flame  speeds  of  Ji  et  al.  The  computed  results 
agree  with  the  experiments  at  Tu=470  K,  and  are  within  the  experimental  uncertainties 
at  Tu=403  K  for  lean  fuel/air  mixtures.  On  the  rich  side,  the  computed  flame  speeds  are 
consistently  higher  than  the  non-linearly  extrapolated  data  from  Ji  et  al .,  while  still  within 
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Figure  43:  Species  profiles  during  the  oxidation  of  n-dodecane/02 /argon  mixtures  in  shock  tubes  as  a 
function  of  time,  at  P=2.25  atm,  0=1.0;  symbols  -  experiments:  Davidson  et  al.  [99];  solid  lines  -  present 
work. 


the  uncertainty  in  the  linearly  extrapolated  flame  speeds,  and  also  agreeing  closely  with  the 
data  of  Kumar  et  al. 
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Figure  44:  OH  concentration  profiles  during  the  oxidation  of  n-dodecane/02/Argon  in  shock  tubes  as  a 
function  of  time  at  P=15  atm,  0=0.5;  symbols  -  experiments:  Vasu  et  al.  [64];  solid  lines  -  present  work. 

Ji  et  al.  [101]  found  that  the  laminar  flame  speeds  of  normal  alkanes  up  to  n-dodecane 
were  similar  at  P=1  atm,  extending  a  previous  such  finding  for  C4  to  C7  alkanes.  Kelley 
et  al.  [102]  measured  laminar  flame  speeds  for  various  n-alkanes  at  higher  pressures,  and 
extended  the  fuel  similarity  to  elevated  pressures.  To  evaluate  the  performance  of  the 
mechanism  at  different  pressures,  and  unburnt  temperatures,  we  seek  to  exploit  the  fuel 
similarity  exhibited  by  n-alkanes. 

The  laminar  flame  speeds  of  n-alkanes/air  mixtures  measured  by  Kelley  et  al.,  at 
Tu=353  K,  and  P=l,  2,  5,  and  10  atm,  and  the  n-dodecane  flame  speeds  computed  using 
the  present  mechanism  are  plotted  in  Figs.  45(b)  and  45(c).  The  computed  results  agree 
closely  with  the  experimental  data  for  all  pressures  examined.  In  addition,  the  agreement 
between  computations  and  experiments  improve  with  increasing  pressure,  a  welcome  trend, 
since  applications  include  engine-like  elevated  pressure  environments. 

The  similarity  exhibited  by  the  normal  alkanes  is  further  explored  by  considering  flame 
speeds  for  n-heptane,  which  is  part  of  the  base  mechanism.  At  P=1  atm,  Ttt=353  K,  and 
0=1.0,  the  computed  flame  speeds  of  n-dodecane/air  and  n-heptane/air  mixtures  are  very 
close,  51.2  cm/s  and  51.6  cm/s,  respectively.  Also,  Fig.  46(a)  shows  very  similar  temperature 
profiles  for  the  two  fuels  as  a  function  of  distance  along  the  flame.  This  is  also  true  of  the 
CO  profiles,  indicative  of  similar  heat  release  for  the  two  fuels,  consistent  with  similar  flame 
speeds.  The  concentrations  of  C2  and  C3  species  from  unimolecular  fuel  decomposition  in 
Fig.  46(b)  show  minor  differences.  However,  upon  examining  the  profiles  of  the  key  radicals 
H,  O,  and  OH  in  Fig.  46(c)  for  the  two  fuels,  it  can  be  concluded  that  these  profiles  are 
virtually  identical,  further  supporting  fuel  similarity  and  consistency  within  the  mechanism. 

n-Dodecane  Summary 

Using  these  validation  studies,  the  base  chemical  mechanism  [59,  80]  has  been  extended  to 
include  the  reaction  pathways  of  n-dodecane  oxidation.  The  resulting  mechanism  has  been 
extensively  validated  for  n-dodecane  oxidation  against  various  experimental  data  sets  that 
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(a)  P=1  atm,  Tu  =  403  K  (b)  P=1  atm,  Tu  =  353  K 


(c)  Higher  pressures,  Tu  =  353  K 

Figure  45:  Laminar  burning  speeds  of  n-dodecane/air  mixtures;  symbols  -  experiments:  (a)  -  Kumar  & 
Sung  [100],  Ji  et  al.  [101],  (b),(c)  -  Kelley  et  al.  [102];  solid  lines  -  present  work. 

included  ignition  delays  and  species  profiles  from  shock  tubes,  and  laminar  burning  velocity 
measurements.  The  observed  similarity  in  ignition  delays  and  laminar  flame  speeds  among 
normal  alkanes  has  been  invoked  to  assess  the  mechanism.  Viewed  collectively,  the  ability 
of  the  present  mechanism  to  predict  different  targets  has  been  found  to  be  satisfactory. 

3.3.3  Extension  To  Methylcyclohexane 

Several  detailed  chemical  mechanisms  have  been  proposed  for  methylcyclohexane  [109,  110, 
111].  Of  specific  interest  here  is  the  chemical  scheme  proposed  by  Pitz  et  al.  [Ill],  combining 
a  low-temperature  mechanism  with  the  high  temperature  mechanism  previously  developed 
by  Orme  et  al.  [110].  This  chemical  scheme  was  validated  against  ignition  delay  time  data 
from  rapid  compression  machine  studies.  The  Pitz  model  has  been  used  as  the  starting 
detailed  mechanism  for  methylcyclohexane  in  the  present  effort,  and  is  reduced  to  a  skeletal 
level  so  that  the  resulting  skeletal  level  mechanism  can  be  combined  with  the  base  model. 
This  is  the  same  base  model  that  was  extended  to  include  aromatics  and  n-dodecane.  The 
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Figure  46:  Analysing  structures  of  n-dodecane  and  n-heptane  premixed  flames  using  the  present  mechanism, 
at  P=1  atm,  Tu  =  353  K,  0=1.0. 

methodology  associated  with  the  extension  to  methycyclohexane  is  discussed  next,  followed 
by  a  demonstration  of  the  performance  of  the  combined  mechanism  for  different  targets. 

Mechanism  Development 

The  reference  mechanism  for  methylcyclohexane  oxidation  from  Pitz  et  al.  [Ill],  which 
will  be  referred  to  as  the  LLNL  mechanism,  has  8807  forward  and  reverse  reactions  and 
999  species.  This  detailed  mechanism  is  first  reduced  to  a  skeletal  level  using  a  multi¬ 
stage  reduction  strategy  put  forth  by  Pepiot-Desjardins  and  Pitsch,  involving  species  and 
reaction  elimination  using  the  DRGEP  approach  [60],  and  chemical  species  lumping  [61]. 
The  database  used  to  carry  out  the  reduction  includes  plug  flow  reactor-like  configurations, 
and  ignition  delays  for  low  to  high  temperatures  (600-1500  K),  pressures  from  1  to  40  atm, 
and  equivalence  ratios  from  0.5  to  1.5. 

The  resulting  skeletal  mechanism  for  methylcyclohexane  consists  of  509  reactions  and 
187  species.  In  Fig.  47,  ignition  delays  computed  using  the  detailed  and  skeletal  schemes 
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Figure  47:  Ignition  delay  times  of  methylcyclohexane/air,  comparing  the  skeletal  level  mechanism  with  the 
detailed  mechanism  results;  solid  lines  -  detailed  mechanism:  Pitz  et  al.  [Ill]  (LLNL). 

are  compared  for  lean  and  stoichiometric  methylcyclohexane/air  mixtures.  The  comparisons 
show  very  good  agreement  for  a  variety  of  temperatures  and  pressures.  A  high  reduction 
ratio  is  achieved  here  (~80%)  and  negligible  error  has  been  introduced  in  the  skeletal  scheme 
relative  to  the  detailed  chemical  mechanism. 

Combination  of  the  skeletal  level  methylcyclohexane  mechanism  and  the  base  mechanism 
is  accomplished  in  a  consistent  fashion,  along  similar  lines  as  that  for  n-dodecane,  described 
in  the  previous  section.  The  ignition  delays  computed  using  the  resulting  mechanism  were 
found  to  be  too  long  compared  to  the  experimental  data  at  all  temperatures  (not  shown 
here).  A  few  changes  have  been  introduced  in  the  combined  mechanism  to  better  describe 
the  kinetics  of  methylcyclohexane  oxidation.  These  changes  will  now  be  discussed. 

Validation  Tests  For  Methylcyclohexane 

The  methylcyclohexane  extension  is  validated  using  ignition  delay  measurements  in  shock 
tubes  and  rapid  compression  machines.  These  experiments  consider  a  wide  range  of  temper¬ 
atures,  pressures,  and  fuel/air  mixture  ratios.  Additionally,  validation  is  performed  using 
OH  and  H2O  profiles  from  shock  tube  measurements,  major  species  measurements  in  Plug 
Flow  Reactors  (PFR),  and  laminar  flame  speed  measurements  at  various  pressures.  The 
complete  list  of  methylcyclohexane  validation  tests,  which  are  discussed  in  further  detail 
below,  appears  in  Table  6. 

Ignition  Delay  Results 

Several  changes  have  been  introduced  in  the  combined  mechanism  to  better  describe  the 
kinetics  of  methylcyclohexane  oxidation,  and  thus  improve  the  agreement  of  the  present 
reaction  scheme  when  compared  to  experimental  ignition  delay  times.  The  reaction  rate 
changes  introduced  are  summarized  here.  First,  bi-radical  pathways  were  removed  from  the 
combined  mechanism  as  suggested  by  Silke  et  al.  [122],  With  this  change,  the  ignition  delays 
at  all  temperatures  showed  better  agreement  with  the  experimental  data.  The  reaction  rate 
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Shock  tubes 

Plug  Flow  Reactor 

Burning  velocity 

Ignition  delays 

Species  profiles 

Vanderover  et  al.  [112] 
Vasu  et  al.  [113] 

Pitz  et  al.  [Ill] 

Hong  et  al.  [114] 

Orme  et  al.  [110] 
Hawthorn  et  al.  [115] 

Vasu  et  al.  [116] 
Hong  et  al.  [114] 

Zeppieri  et  al.  [117] 

Kumar  et  al.  [118] 
Singh  et  al.  [119] 

Ji  et  al.  [120] 

Kelley  et  al.  [121] 

Table  6:  Validation  cases  for  methylcyclohexane  undertaken  in  the  present  study. 


for  the  decomposition  reaction: 


C5H9  ->  C2H3  +  C3H6 

was  found  to  be  sensitive  for  high  temperature  ignition,  and  this  rate  was  adopted  from 
Tsang  et  al.,  leading  to  faster  ignition  at  high  temperatures. 

Improvements  in  NTC  region  ignition  delays  were  achieved  by  updating  reaction  rates 
for  the  following  key  reactions  based  on  recent  theoretical  studies.  The  ignition  delays  in  the 
NTC  region  were  found  to  be  sensitive  to  (i)  ring  opening  of  methylcyclohexyl  radicals,  for 
instance  as  depicted  in  Fig.  48(a),  (ii)  reactions  involving  isomerization  of  peroxy  radicals  to 
hydroperoxy  radicals,  as  in  Fig.  48(b),  and  (iii)  reactions  for  formation  of  epoxides  (cyclic 
ether)  and  OH  from  the  hydroperoxy  radical  as  shown  in  Fig.  48(c).  Recent  theoretical 
studies  by  Sirjean  et  al.  have  focussed  on  similar  reactions  for  cyclohexyl  radicals  [123],  cy¬ 
clohexyl  peroxy,  and  hydroperoxy  radicals  [124].  In  the  present  chemical  model,  the  reaction 
rates  for  the  methylcyclohexyl  ring  opening  and  closure  are  adapted  from  those  of  cyclo¬ 
hexyl  radicals  [123],  Similarly,  the  rate  for  the  formation  of  epoxide  from  methylcyclohexyl 
hydroperoxy  radical,  and  the  isomerization  of  methylcyclohexyl  peroxy  radical  is  adapted 
from  that  of  the  cyclohexyl  counterparts  [124].  The  reverse  isomerization  rate  is  computed 
from  thermodynamic  properties.  These  rate  changes  help  achieve  better  agreement  with 
experimental  data  in  the  NTC  region. 

Methylcyclohexane/ Air  Ignition  Delays 

The  performance  of  the  mechanism  for  ignition  delays  in  methylcyclohexane/air  mixtures 
is  demonstrated  by  comparing  against  experimental  ignition  delay  time  data  from  the  shock 
tube  experiments  of  Vasu  et  al.  [113]  at  (f>  =  1.0  and  P=20  and  45  atm.  Comparisons  are 
also  performed  with  the  data  of  Vanderover  and  Oehlschlaeger  [112]  at  <f>  =  0.25,0.5,1.0 
and  P=12  and  50  atm.  Very  good  agreement  is  obtained  for  ignition  delay  times  between 
the  present  scheme  and  the  experimental  data  in  Fig.  49  for  low  through  high  temperatures 
as  a  result  of  the  changes  introduced  above.  Also  shown  in  these  figures  are  comparisons 
against  ignition  delay  times  measured  in  a  Rapid  Compression  Machine  by  Pitz  et  al.  [Ill] 
at  (j)  =  1.0,  with  P  equal  to  either  15  or  20  atm.  The  present  scheme  captures  the  pressure 
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(c)  Hydroperoxy  radical  forming  cyclic  ether  +  OH 
Figure  48:  Methylcyclohexane  reactions  that  are  important  in  the  NTC  ignition  regime. 

dependence  of  the  ignition  delays  as  well  the  variation  with  fuel  mixture  ratios  with  good 
accuracy. 

Methylcyclohexane/ 02/Ar  Ignition  Delays 

Ignition  delay  time  measurements  have  been  made  for  methylcyclohexane/02/Ar  mix¬ 
tures  at  near-atmospheric  pressures  by  Hong  et  al.  [114]  for  cj)  =  0.5  and  cf)  =  1.0  mixtures. 
A  comparison  of  the  predicted  ignition  delay  times  and  the  experimental  data  is  shown  in 
Figs.  50(a)  and  50(b).  The  computed  ignition  delays  show  good  agreement  with  the  ex¬ 
perimental  data  considered  in  Fig.  50(a),  but  shorter  ignition  delays  than  experiments  for 
T  <  1200  K  in  Fig.  50(b). 

Other  ignition  delay  time  measurements  for  methylcyclohexane/02/Ar  mixtures  at  sim¬ 
ilar  pressures  come  from  Vasu  et  al.  [113],  Orrne  et  al.  [110],  and  Hawthorn  and  Nixon  [115]. 
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Figure  49:  Ignition  delay  times  of  methylcyclohexane/air  mixtures:  Symbols  -  experimental  data  from 
Vanderover  and  Oehlschlaeger  [112],  Pitz  et  al.  [Ill],  Vasu  et  al.  [113];  Solid  lines  -  ignition  delays  computed 
using  the  present  reaction  scheme. 
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Figure  50:  Ignition  delay  times  of  methylcyclohexane/02/Ar  mixtures:  Symbols  -  experimental  data  from 
(a,b,c)  Hong  et  al.  [114],  (c)  Vasu  et  al.  [113],  Orme  et  al.  [110],  Hawthorn  and  Nixon  [115];  Solid  lines  - 
ignition  delays  computed  using  the  present  reaction  scheme. 


In  Fig.  50(c),  these  data,  along  with  those  of  Hong  et  al.  are  scaled  to  a  common  fuel  concen¬ 
tration  Xmch  =  0.381%,  Xo2=2%,  balance  Ar  (</>  =  1.0),  P=1.5  atm,  using  the  correlation 
proposed  by  Vasu  et  al.  [113], 

~  ~  p— 0.98±0.13v— 0.82±0.08  j1.47±0.09 

'ign  tx.  Jr  ^MCH  P 

The  ignition  delays  computed  at  these  conditions  using  the  present  scheme  fall  within  the 
variability  in  the  scaled  experimental  data  in  Fig.  50(c). 

Results  For  OH/H2O  Time- Histories 

Vasu  et  al.  [116]  measured  OH  profiles  during  methylcyclohexane/02/Ar  oxidation,  with 
Xmch  =  1000  ppm,  Xo2  =  0.021,  balance  Ar,  at  an  equivalence  ratio  of  cf>  =  0.5  and  a 
pressure  of  P=15  atm.  The  OH  profiles  predicted  using  the  present  reaction  scheme  are 
compared  with  the  experimentally  measured  profiles  of  Vasu  et  al.  in  Fig.  51.  The  agreement 
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Figure  51:  Methylcyclohexane/02/Ar  oxidation:  (a-c)  OH  profiles  at  Xmch  =  1000  ppm,  Xo2  =  0.021, 
balance  Ar,  ((f)  =  0.5)  P=16  atm  (d)  Ignition  delay  times  corresponding  to  (i)  the  experimental  conditions 
of  (a-c),  and  (ii)  Xmch  =  750  ppm,  (f>  =  0.5,  P=16  atm;  Symbols  -  experimental  data  from  Vasu  et  al.  [116]; 
Solid  lines  -  results  computed  using  the  present  reaction  scheme. 


is  satisfactory  across  the  conditions  being  considered,  with  the  peak  OH  value  and  the  rise 
of  OH  profiles  predicted  reasonably  well.  The  ignition  delays  were  also  measured  by  Vasu  et 
al.  at  the  same  conditions  at  which  the  OH  profiles  were  obtained,  along  with  another  set  of 
ignition  delay  data  at  Xmch  =  750  ppm,  4>  =  0.5,  and  P=15  atm.  These  data  are  compared 
to  the  ignition  delay  times  computed  using  the  present  chemical  model  in  Fig.  51(d),  and 
show  good  agreement. 

Hong  et  al.  [114]  also  measured  OH  and  H2O  profiles  during  methylcyclohexane/Oo/Ar 
oxidation,  at  </>  ~  0.85  and  P=2.2  atm  in  a  shock  tube  facility.  The  computed  species  profiles 
are  compared  with  the  experimentally  measured  profiles  of  Hong  et  al.  in  Fig.  52.  While 
the  predictions  show  good  agreement  with  the  experimental  peak  values  and  rise  of  OH  and 
H2O  time-histories,  at  T=1359  K  the  predicted  rise  of  profiles  occurs  earlier  than  in  the 
experiments.  This  early  rise  of  the  OH/H2O  profiles  was  also  observed  when  comparing  the 
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Figure  52:  Species  time-histories  during  methylcyclohexane/02/Ar  oxidation.  Conditions:  (1)  T  =  1359  K, 
P=2.2  atm,  380  ppm  MCH/4200  ppm  02/Ar,  <j>  =  0.95,  (2)  T  =  1453  K,  P=2.2  atm,  340  ppm  MCH/4200 
ppm  02/Ar,  cj)  =  0.85,  (3)  T  =  1541  K,  P=2.1  atm,  320  ppm  MCH/4200  ppm  02/Ar,  ij>  =  0.80;  Symbols  - 
experimental  data  from  Hong  et  al.  [114];  Dashed  lines  -  OH  and  H20  profiles  computed  using  the  present 
reaction  scheme. 


jetSurF  mechanism  predictions  with  these  experimental  data  (Fig.  15  in  Hong  et  al.  [114]), 
and  this  warrants  further  investigation. 

Laminar  Flame  Speed  Results 

Several  experimental  studies  have  measured  laminar  flame  speeds  of  methylcyclohex- 
ane/air  mixtures  using  different  measurement  techniques.  At  atmospheric  pressure  and  an 
unburnt  temperature  of  Tu=353  K,  Ji  et  al.  [120]  determined  flame  speeds  in  a  counterflow 
configuration.  Kelley  et  al.  [121]  measured  propagation  speeds  of  spherically  expanding 
methylcyclohexane/air  flames  at  several  pressures  with  Tu=353  K.  Similar  techniques  were 
used  by  Singh  et  al.  [119]  to  measure  flame  speeds  at  P=1  atm  and  Tu=400  K.  Also,  Ku¬ 
mar  &  Sung  [118]  used  a  counterflow  twin- flame  technique  to  determine  flame  speeds  of 
methylcyclohexane/air  mixtures  at  these  conditions. 

The  flame  speeds  computed  using  the  present  reaction  scheme  are  compared  with  these 
experimental  data  in  Fig.  53.  At  atmospheric  pressure,  the  flame  speeds  at  Tu=400  K  in 
Fig.  53(a)  agree  closely  with  the  data  of  Singh  et  al,  while  at  Tu  =  353  K  in  Fig.  53(b),  the 
computed  curves  agree  closely  with  the  data  of  Kelley  et  al.  The  degree  of  agreement  with 
the  Kelley  et  al.  data  remains  consistently  good  at  the  higher  pressures  as  well.  In  all,  the 
ability  of  the  reaction  model  to  predict  laminar  flame  speeds  is  demonstrated  to  be  good, 
and  within  the  variability  of  the  experimental  flame  speed  data  determined  using  different 
measurement  techniques. 

Species  Profiles  In  Plug  Flow  Reactors 

High  temperature  oxidation  of  methylcyclohexane/air  mixtures  was  studied  by  Zeppieri  et 
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Figure  53:  Laminar  flame  propagation  speeds  in  methylcyclohexane/air  mixtures:  Symbols  -  experimental 
data  from  Kumar  &  Sung  [118],  Deepti  Singh  &  Qiao  [119],  Ji  et  al.  [120],  Kelley  et  al.  [121];  Solid  lines  - 
flame  speeds  computed  using  the  present  reaction  scheme. 


al.  [117]  in  the  Princeton  Turbulent  Flow  Reactor.  They  measured  concentrations  of  methyl- 
cyclohexane  (MCH),  as  well  as  typical  major  and  minor  species  formed  during  oxidation. 
The  experimental  data  was  obtained  at  P=1  atm,  T  =  1160  K,  initial  MCH  concentration 
of  1815  ppm,  and  fuel/air  mixture  ratio  of  <j>  =  1.3.  The  experimental  configuration  is 
modeled  as  a  isobaric  reactor  in  the  present  work. 

A  comparison  between  the  species  concentration  time  evolution  predicted  using  the 
present  scheme  and  the  experimental  data  is  presented  in  Fig.  54.  The  change  in  temper¬ 
ature  as  a  function  of  time  is  accurately  represented  by  the  present  scheme  in  Fig.  54(b). 
The  fuel  decay  is  faster  than  the  experiments  in  Fig.  54(a),  which  is  also  consistent  with  the 
discrepancies  between  the  predictions  in  ethene  and  ethane  profiles  and  the  experiments. 
An  improvement  in  the  fuel  decay  profile  is  therefore  expected  to  lead  to  better  agreement 
in  C2H4  and  C2H5  profiles  as  well.  The  evolution  of  CO2  has  been  predicted  accurately  in 
Fig.  54(d).  However,  some  discrepancies  are  seen  in  CO,  C3H6  and  CH4  concentrations  in 
Figs.  54(c),  54(g),  and  54(h). 

Methylcyclohexane  Summary 

With  the  extensive  validation  tests  presented  here,  it  can  be  concluded  that  the  present 
methylcyclohexane  oxidation  model  is  able  to  reproduce  experimental  ignition  delay  times, 
OH/H2O  profiles,  major  species,  and  laminar  flame  speeds  with  reasonable  accuracy. 
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Figure  54:  Major  species  time-evolution  during  methylcyclohexane/air  oxidation  in  a  Plug  Flow  Reactor  at 
P=1  atm,  T  =  1160  K,  X^ch  =  1815  ppm,  (f>  =  1.3:  Symbols  -  experimental  data  from  Zeppieri  et  al.  [117]; 
Solid  lines  -  species  profiles  computed  using  the  present  reaction  scheme. 
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Re- Validation  Of  Individual  Components  In  The  Surrogate  Mechanism 

With  the  addition  of  the  reaction  pathways  of  methylcyclohexane,  several  aromatics  (es¬ 
pecially  m-xylene),  and  n-dodecane  to  the  base  model,  the  combined  chemical  mechanism 
can  be  used  to  describe  the  oxidation  of  the  jet  fuel  surrogate  proposed  in  Table  2.  The  final 
step  of  this  process  consists  of  re- validating  the  different  surrogate  fuel  components  now  that 
they  have  been  combined  to  form  the  final  surrogate  mechanism.  This  re-validation  is  used 
to  evaluate  any  non-linear  changes  that  might  have  occurred  because  of  the  added  reaction 
modules.  The  results  for  the  different  test  cases  for  the  individual  surrogate  components 
were  recomputed  using  this  reaction  scheme,  and  little  changes  were  observed  from  the 
previous  results  (not  shown  here).  Most  of  the  differences  resulted  in  improved  agreement 
with  the  relevant  experimental  data. 

3.4  Performance  Of  The  Model  Jet  Fuel  Surrogate 

The  capabilities  of  the  model  JP-8  surrogate  proposed  in  Table  2  are  now  evaluated.  Global 
kinetic  targets  such  as  ignition  delay  times  and  laminar  flame  speeds  computed  using  the 
present  reaction  scheme  are  compared  with  experimental  data  for  real  JP-8  fuels  in  this 
section. 

Ignition  Delay  Results 

Figure  55  shows  a  comparison  of  ignition  delays  predicted  using  the  present  scheme  with 
experimental  measurements  from  Vasu  et  al.  [125].  The  experiments  were  conducted  in  a 
shock  tube  at  P=20  atm,  using  fuel/air  mixture  ratios  of  (p  =  0.5  and  1.0.  Agreement  of  the 
model  and  the  experiments  is  excellent  for  both  of  the  equivalence  ratios  that  are  considered 
here. 

Laminar  Flame  Speed  Results 

Ji  et  al.  [126]  measured  laminar  flame  speeds  of  JP-7  and  JP-8  fuels  at  atmospheric 
pressure  and  at  an  unburnt  temperature  of  Tu  =  403  K.  The  computed  flame  speeds  for 
the  model  surrogate  are  compared  with  this  experimental  data  in  Fig.  56.  The  flame  speed 
predictions  closely  follow  the  JP-8  flame  speeds  at  all  fuel/air  ratios.  It  can  also  be  seen 
that  the  experimentally  determined  flame  speeds  of  n-dodecane  are  approximately  equal  to 
those  of  the  jet  fuels,  justifying  the  use  of  this  component  as  the  significant  alkane  within 
the  model  surrogate. 

The  ability  of  the  present  scheme  to  predict  global  oxidation  characteristics  of  real  jet 
fuels,  such  as  ignition  delays  and  laminar  flame  speeds,  can  be  concluded  to  be  excellent. 
The  chemical  model  will  be  further  assessed  in  other  configurations  of  interest  in  future 
work. 
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Figure  55:  Ignition  delay  times:  Symbols  -  experimental  data  from  Vasu  et  al.  [125];  Solid  lines  -  predictions 
using  the  present  reaction  scheme  with  the  model  surrogate  proposed  in  Table  2  as  the  fuel. 
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Figure  56:  Laminar  flame  speeds:  Symbols  -  experimental  data  from  Ji  et  al.  [126],  Solid  lines  -  flame  speed 
predictions  using  the  present  reaction  scheme  with  the  model  surrogate  proposed  in  Table  2  as  the  fuel. 
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4  Conclusions 


This  project  has  addressed  physical  modeling  needs  in  the  areas  of  both  sub- filter  LES 
model  development  and  chemical  mechanism  formulation  and  validation.  The  resulting 
models  are  expected  to  directly  influence  and  improve  the  next  generation  of  high  fidelity 
reactive  flow  simulations.  Additionally,  the  two  project  components  naturally  interface 
within  a  comprehensive  flow  modeling  framework.  For  example,  chemical  mechanisms  can 
be  specifically  tailored  for  individual  modeling  studies  using  the  surrogate  framework.  A 
highly  detailed  n-dodecane  mechanism  might  be  reduced  and  combined  with  a  particular 
aromatic  mechanism  for  use  in  a  DNS  study  of  soot  that  relies  on  finite  rate  chemistry. 
Conversely,  the  same  mechanism  might  be  left  in  an  unreduced  form  and  combined  with 
several  aromatic  mechanisms  for  use  in  an  LES  based  on  the  flamelet  approach.  The 
surrogate  framework  enables  this  level  of  tailoring  of  the  chemistry,  so  that  the  modeling 
approach  can  be  matched  to  the  application. 

Progress  in  the  area  of  sub-filter  LES  modeling  similarly  sets  the  stage  for  combustion 
model  adaptation.  The  individual  advances  in  non-premixed  and  premixed  combustion 
modeling  might  be  applied  directly  in  simulations  of  engines  characterized  by  a  single  com¬ 
bustion  regime.  Conversely,  in  more  complex  settings  where  multiple  combustion  regimes 
are  expected,  the  multi-regime  flamelet  approach  can  be  employed.  The  advantages  of 
the  new  single  regime  developments  can  then  be  leveraged  where  needed,  using  a  chemical 
mechanism  that  is  tailored  to  describe  either  flame  extinction  or  burning  velocities. 

In  future  work  these  newly  developed  models  will  continue  to  be  advanced  and  integrated 
into  Stanford’s  computational  infrastructure.  Furthermore,  they  will  made  available  to  the 
Air  Force  Office  of  Scientific  Research  upon  request.  The  next  step  in  the  process  of  fully 
leveraging  their  potential  is  to  integrate  their  capabilities  in  the  simulation  of  complex 
geometry  engines.  This  is  the  planned  topic  of  future  studies. 
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